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1.0  INTRODUCTION 

Electromagnetic  coupling  and  shielding  problems  have  traditionally 
been  difficult  to  treat  with  analytical  or  numerical  methods  because  of 
the  failure  of  these  methods  to  adequately  resolve  the  effects  of  shield 
apertures,  curvatures,  corners,  and  internal  contents.  Usually,  only 
relatively  simple  geometries  of  shields  and  shield  openings  are  studied 
in  an  attempt  to  gain  insight  into  the  key  coupling  mechanisms,  and  to 
allow  a rough  estimate  of  the  coupling  for  more  complicated  and  realistic 
problems.  A method  for  the  direct  modeling  and  solution  of  realistic 
problems  would  eliminate  the  need  for  intuition  in  applying  simple  models, 
and  would  greatly  increase  the  accuracy  of  the  ultimate  result. 

This  research  program  investigated  the  application  of  a new  approach 
for  the  direct  modeling  of  electromagnetic  interaction  problems:  the 
finite-difference,  time-domain  (FD-TD)  solution  of  Maxwell's  equations. 

The  FD-TD  method  treats  the  irradiation  of  a structure  as  an  initial  value 
problem.  At  t = 0,  a plane  wave  source  of  frequency,  f,  is  assumed  to  be 
turned  on.  The  propagation  of  waves  from  this  source  is  simulated  by 
solving  a finite-difference  analog  of  the  time-dependent  Maxwell's  equations 
on  a lattice  of  cells,  including  the  structure.  Time-stepping  is  continued 
until  the  sinusoidal  steady  state  is  achieved  at  each  cell.  The  field 
envelope,  or  maximum  absolute  value,  during  the  final  half  wave-cycle  of 
time-stepping  is  taken  as  the  magnitude  of  the  phasor  of  the  steady-state 
field. 

This  method  has  two  key  advantages  relative  to  available  modeling 
approaches.  First,  it  is  simple  to  implement  for  complicated  metal/ 
dielectric  structures  because  arbitrary  electrical  parameters  can  be 
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assigned  to  each  lattice  cell  using  a data  card  deck.  Second,  its  computer 
memory  and  running  time  requirement  is  not  prohibitive  tor  many  complex 
structures  of  interest. 

This  report  first  reviews  available  numerical  techniques  for  the 
solution  of  electromagnetic  coupling/shielding  problems.  Then,  the  basic 
elements  of  the  FD-TD  method  are  introduced,  with  detailed  derivations 
where  appropriate,  and  with  examples  of  prior  computed  results  for 
dielectric  structures  to  establish  the  expected  level  of  accuracy  for 
general  structures.  The  last  section  of  this  report  will  detail  the  com- 
puted results  for  the  metal  geometries  modeled  in  this  program  effort. 

Full  listings  of  computer  programs  employed  in  this  research  effort  are 
provided  in  the  Appendix. 
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2 . 0 REVIEW  OF  AVAILABLE  NUMERIC^.  TECHNIQUES 


The  coupling  of  electromagnetic  fields  to  the  interior  of  an  arbitrary 
conducting  or  dielectric  structure  has  been  approached  using  both  frequency 
domain  and  time  domain  numerical  techniques.  This  section  reviews  pub- 
lished work  in  this  area  and  discusses  the  problems  inherent  in  present 
approaches. 

2.  1 Frequency  Domain  Techniques 

Frequency  domain  methods  are  based  upon  the  assumption  of  an  exp(j2TTft) 
time  dependence  in  the  fundamental  Maxwell's  equations.  In  general,  methods 
of  this  type  derive  a set  of  linear  equations  for  either  field  variables  or 
field  expansion  coefficients,  and  then  solve  the  linear  system  with  a suit- 
able matrix  inversion  scheme. 

Almost  all  frequency  domain  techniques  can  be  placed  in  the  following 
three  classes. 

1.  Electromagnetic  field  expansions,  in  terms  of  either 

a.  Normal  modes  of  the  structure,^ 

2 

b.  Analytic  continuation  of  free  space  modes, 

c.  Normal  modes  of  the  structure  matched  to  aperture  fields 

3 

determined  using  a quasi-static  approximation, 

d.  Normal  modes  of  the  structure  matched  to  aperture  fields 
resulting  from  an  exterior  region  expansion  of  free 

4 

space  modes. 

Using  these  techniques,  the  solution  is  achieved  by  enforcing  the  boundary 
conditions  for  the  fields  at  a sufficient  number  of  points  to  specify  the 
surfaces  of  the  structure,  and  to  obtain  a set  of  simultaneous  equations 
for  the  modal  coefficients. 

2.  Integral  equation  solutions,  set  up  at  either 

5 

a.  Surfaces  of  the  structure, 

b.  Internal  volumes  of  the  structure,^ 

c.  Surfaces  of  a structure  with  aperture.^ 
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Here,  the  solution 'is  achieved  by  enforcing  the  integral  eguation  at  a 
sufficient  number  of  points  to  specify  the  surfaces  or  volume  of  the  struc- 

Q 

ture,  using  the  method  of  moments  to  obtain  a set  of  simultaneous  eguations 
for  the  electromagnetic  fields  at  eaci)  enforcement  point. 

3.  Variational  solutions,  employing  specifically  the  finite  element 

Q 

solution  of  the  Helmholtz  eguation  within  an  unbounded  region.  Here,  the 
solution  is  obtained  by  enforcing  the  variational  principal  to  obtain  a set 
of  simultaneous  equations  for  either  field  variables  or  expansion  coeffi- 
cients . 

In  principle,  the  accuracy  of  frequency  domain  methods  is  excellent  if 
a sufficiently  large  set  of  simultaneous  equations  is  solved.  However,  ea-^h 
method  may  have  one  of  two  problems  when  complicated  structures  are  con- 
sidered. First,  excessive  computer  storage  may  be  required.  Second,  an 
excessively  complex  derivation  of  the  matrix  elements  may  be  required.  We 
now  consider  these  possibilities  in  the  context  of  modeling  the  interior 
fields  of  an  arbi trary  structure  having  dimensions  on  the  order  of  one 
wavelength. 

To  illustrate  the  problem  of  excessive  computer  storage  requirements, 

consider  the  integral  equation  technique  of  Ref.  6.  Using  the  authors' 

maximum  allowed  solution  point  spacing  of  1/4  wavelength,  a total  of  about 
3 

3x4  = 192  equations  must  be  solved  to  determine  the  3 rectangular  com- 
ponents of  the  electric  field  (or  magnetic  field)  at  each  of  the  approxi- 

3 2 

mately  4 points  of  the  structur-e.  For  this  case,  a maximum  of  192  == 

37,000  interaction  coefficients  must  be  stored  in  the  computer  to  allow 

inversion  of  the  system  matrix,  assuming  no  matrix  bandl imi ting.  Yet,  a 

1/4  wavelength  resolution  might  be  much  too  coarse  to  resolve  needed  details 

of  the  structure.  If  a solution  point  spacing  of  1/8  wavelength  is  selected, 

3 

the  number,  N,  of  simultaneous  equations  increases  to  3 x 8 = 1536,  and  a 
7 fi 

total  of  N 2.3  X 10  coefficients  must  be  stored,  a 64-fold  increase  over 
the  previous  case.  It  is  seen  that,  for  a fixed  structure  geometry,  the 
required  computer  storage  varies  approximately  as  the  inverse  of  the  sixth 
power  of  the  resolution.  Conversely,  for  a fixed  resolution,  the  required 
storage  varies  as  the  sixth  power  of  the  structure's  characteristic  dimen- 
sion. This  extremely  rapid  rate  of  increase  limits  the  application  of 


many  frequency  domain  techniques  to  either  two-dimensional  irradiation  prob- 
lems (infinitely  long  cylinders  of  constant  cross  section),  or  simple  three- 
dimensional  problems  (bodies  of  revolution,  wire  grid  models  of  conducting 
surfaces),  which  can  be  solved  with  relatively  few  simultaneous  equations. 

To  illustrate  the  problem  of  an  excessively  complex  derivation  of  the 
matrix  elements,  we  consider  the  technique  of  Ref.  2.  For  each  dielectric 
medium  of  the  arbitrary  structure,  we  would  have  to  assume  two  field  expan- 
sions: one  inside  the  medium,  and  one  outside  the  medium.  Field  matching 
would  be  performed  at  enough  points  along  the  surface  of  each  medium  so  that 
its  shape  would  be  outlined.  At  each  field  matching  point,  an  analytic  con- 
tinuation of  the  interior  and  exterior  fields  would  be  formulated.  Continua- 
tion of  the  interior  fields  for  some  media  would  require  multiple  individual 
continuations,  to  account  for  elongated  geometries,  implying  multiple  sum- 
mations of  spherical  Bessel  functions.  Any  change  of  the  structure  geometry 
would  require  a recalculation  of  virtually  every  matrix  element  involving 
new  analytic  continuations.  Although  the  size  of  the  matrix  obtained  with 
this  method  is  not  excessive,  the  derivation  of  the  matrix  elements  is  com- 
plicated, with  the  complexity  of  the  derivation  a function  of  the  inhomo- 
geneities and  geometry  of  the  structure.  This  implies  that  a long  program 
development  time  is  required  for  each  new  problem. 

2.2  Previous  Time  Domain  Techniques 

Most  previous  time  domain  techniques  can  be  placed  in  the  following 
four  classes : 

1.  Inverse  Fourier  transforms  of  frequency  domain  solutions,  either 

a.  Obtained  from  previously  calculated  frequency  domain  solutions^ 
possibly  via  the  fast  Fourier  transform, 

b.  Derived  in  the  form  of  a time  convolution.^^ 

2.  Transformed  s-plane  methods,  for  example  using 

a.  Field  expansions,  matching  at  apertures,  and  the  method  of 

12 

moments , 

1 3 

b.  Singularity  expansion  methods. 
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3.  Integral  equation  solutions,  for  example,  of  the  magnetic  field, 
yielding  the  surface  current  density  induced  on  bodies  of  revolution 

4.  Transmission  line  models^^  and  antenna  models^^  of  induced  currents 
and  coupling  at  apertures. 

The  first  class  of  time  domain  techniques  is  seen  to  suffer  the  com- 
puter storage  problems  inherent  in  frequency  domain  methods,  with  the  added 
complexity  of  requiring  Fourier  transform  computational  processing.  For  the 
second  class,  field  expansion  and  solution  via  the  method  of  moments  leads 
again  to  computer  storage  problems,  while  the  singularity  expansion  method 
is  useful  for  determining  the  exterior  surface  current  and  scattered  field, 
rather  than  the  internal  field  distribution.  The  third  class  of  techniques 
requires  computer  storage  of  the  surface  current  and  its  time  derivative  at 
all  surface  points,  for  all  times  between  the  start  of  irradiation  and  the 
observation  time.  This  again  leads  to  computer  storage  problems,  as  did 
method  of  moments  techniques  discussed  before.  The  last  group  of  techniques 
is  suitable  principally  for  thin  diameter,  cylinder-like  structures,  and 
cannot  account  for  details  of  the  interiors  of  the  structures. 

In  summary,  these  four  classes  of  time  domain  techniques  are  probably 
unsuited  for  application  to  shielding  and  coupling  problems  requiring  a good 
amount  of  detail  of  the  interior  of  the  structure  of  interest.  Computer 
storage  is  again  a significant  problem,  as  it  was  for  the  available  fre- 
quency domain  techniques.  The  FD-TD  method,  discussed  next,  still  has  the 
important  advantage  of  allowing  excellent  resolution  of  the  details  of  the 
interior  of  a structure  without  exhausting  the  storage  available  in  large, 
modern-day  computers. 
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3.0  THE  FD-TD  METHOD 


This  section  will  discuss  the  following  topics  relevant  to  the  use  of 
the  FD-TD  method:  1)  a brief  history  of  the  development  of  the  method; 

2)  the  basic  ideas  behind  the  method;  3)  its  advantages  over  the  present 
commonly  used  approaches;  4)  computational  details;  and  5)  past  usages  and 
results . 

3 . 1 Development  of  the  Method 

The  basic  finite-difference  lattice  structure  and  time  stepping  algor- 
ithm for  the  FD-TD  method  was  presented  by  Yee  in  a 1966  paper. Yee, 
however,  was  unable  to  solve  the  problems  of  reflection  of  outgoing  waves 
at  the  lattice  truncations  and  the  generation  of  a long  pulse  or  continuous 
wave  without  elongating  the  lattice.  Therefore,  his  results  were  only  of 
limited  use  for  practical  coupling  and  shielding  problems. 

The  principal  investigator  of  the  present  contract  performed  the  basic 

research  needed  to  fully  develop  the  FD-TD  method  as  part  of  his  doctoral 

dissertation  at  Northwestern  University  (1975).  The  developments  included 

satisfactory  approximations  to  the  free  space  condition  at  the  lattice 

truncations,  and  the  simulation  of  a long  duration  pulse  or  continuous  wave 

incident  on  tiie  structure  of  interest.  The  basic  concepts  of  the  FD-TD 

method  and  the  initial  description  of  the  practical,  wave  scattering  algor- 

1 9 

ithm  were  presented  in  a paper  in  August  1975.  This  paper  applied  the 
method  to  solve  for  the  standing  wave  pattern  within  circular  dielectric 
cylinders  subjected  to  microwave  irradiation.  Additional  problems  of  this 
type,  dealing  with  a square  dielectric  cylinder  and  a dielectric  sphere  were 
solved  using  the  FD-TD  method  as  part  of  the  dissertation.  In  each  of  these 
cases,  the  availability  of  exact  analytic  solutions  permitted  comparison 
with  the  results  of  the  FD-TD  method,  and  the  determine" ion  of  an  error 
bound. 

Recognizing  that  the  FD-TD  method  is  ideal  for  modeling  dielectric 
structures  of  considerable  complexity  and  inhaiiogeneity , the  principal 
investigator  applied  the  method  to  an  important  and  previously  incalculable 
problem,  the  solution  of  the  electromagnetic  fields  and  induced  heating 
potential  within  the  human  eye  subjected  to  microwave  irradiation.  The 
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FD-TD  method  permitted  the  modeling  of  the  human  eye  and  the  various  tissues 

of  the  bony  orbit  surrounding  the  eye  in  great  detail.  The  results  of  this 

20 

study  were  published  in  a paper  in  November  1975.  Experimental  work  per- 
formed by  others  using  rabbits  has  tended  to  confirm  the  predictions  of  the 
FD-TD  eye  model . 

3.2  Basic  Ideas 

3.2.1  Wave  Tracking 

The  FD-TD  method  achieves  its  flexibility  by  programming  the  two  funda- 
mental equations  of  classical  electromagnetic  theory;  Maxwell's  time  depen- 
dent curl  equations.  Using  this  method,  we  can  model  the  propagation  of  an 
arbitrary  pulsed  or  continuous  wave  EM  field  into  a volume  of  space  con- 
taining either  a dielectric  or  conducting  structure.  By  time-stepping,  i.e., 
repeatedly  solving  the  finite-difference  analog  of  the  curl  equations  at 
each  point  of  a space  lattice  containing  the  structure  of  interest,  we 
actually  track  the  incident  wave  as  it  first  propagates  to  the  structure, 
and  then  interacts  with  it  in  some  way  (surface  current  excitation,  diffus- 
ion, penetration,  etc,).  Wave  tracking  is  completed  for  pulsed  irradiation 
when  the  desired  early  or  late  time  behavior  is  observed;  for  sinusoidal 
irradiation,  the  end  point  is  the  attainment  of  the  sinusoidal  steady  state. 

Time-stepping  for  the  FD-TD  method  is  accomplished  by  what  is  termed  an 
explicit  finite-difference  procedure.  Here,  the  value  of  an  electromagnetic 
field  component  at  the  latest  time  step  is  computed  using  only  field  quan- 
tities found  during  the  previous  time  step,  and  stored  in  the  computer 
memory.  For  example,  a particular  electric  field  component,  E^,  to  be 
evaluated  at  point,  P,  of  the  finite-difference  lattice,  is  computed  using 
the  stored  value  of  E^  at  P and  the  stored  values  of  the  magnetic  field  com- 
ponents, and  H^,  at  the  lattice  points  immediately  adjacent  to  P.  In  this 
way,  the  unknown  field  quantity  is  an  explicit  function  of  known  field  quan- 
tities. Thus,  no  simultaneous  equations  are  needed  to  compute  the  fields  at 
the  latest  time  step.  Further,  computation  can  proceed  one  lattice  point  at 
a time,  and  the  new  field  value  at  each  point  can  be  placed  immediately  in 
memory. 
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3.2.2  Modeling  the  Structure  of  Interest 

The  finite-difference  formulation  of  the  FD-TD  method  allows  the  simple 
and  straightforward  modeling  of  arbitrary  dielectric/conducting  structures. 
This  is  because  the  space  containing  the  model  of  the  structure  is  divided 
up  into  discrete  volumes,  or  unit  cells.  The  simplest  case  is  that  of  a 
cubic  unit  cell,  which  results  in  a cubic  lattice  approximation  of  the 
geometry.  For  this  case,  the  structure  of  interest  is  mapped  into  the  space 
lattice  by  first  choosing  the  dimensions  of  the  unit  cell,  and  then  assign- 
ing appropriate  values  of  electrical  permittivity  and  conductivity  to  each 
unit  cell  of  the  lattice.  Thus,  inhomogeneities  or  fine  details  of  the 
structure  can  be  modeled  with  a maximum  resolution  of  one  unit  cell;  thin 
surfaces  can  be  modeled  as  infinitely  thin,  stepped-edge  sheets.  The  com- 
puter program  is  written  so  that  only  a set  of  data  cards  is  required  to 
specify  the  complete  geometry  and  dielectric  characteristics  of  an  arbitrary 
structure.  No  special  handling  of  electromagnetic  boundary  conditions  at 
media  interfaces  is  required  because  the  Maxwell  curl  equations  generate 
these  conditions  in  a natural  way  by  themselves.  Therefore,  the  basic  com- 
puter program  need  not  be  modified  to  change  from  structure  to  structure, 
assuming  that  the  lattice  volume  used  is  sufficient  to  fully  contain  each 
structure. 

3.2.3  The  Lattice  Truncation  Conditions 

A basic  problem  with  any  finite-difference  solution  of  Maxwell's  equa- 
tions in  an  unbounded  region  is  the  treatment  of  the  field  vector  components 
at  the  lattice  truncation.  Because  of  computer  storage  limitations,  the 
lattice  must  terminate  close  to  the  model  structure  in  a region  where  the 
nature  of  the  scattered  wave  is  not  clearly  known.  Proper  truncation  of  the 
lattice  requires  that  any  outgoing  wave  disappears  at  the  lattice  boundary 
without  reflection  during  the  continuous  time  stepping  of  the  algorithm. 
Improper  truncation  would  cause  error  for  all  time  steps  after  the  spurious 
wave  reflections  return  to  the  vicinity  of  the  model  structure. 

The  FD-TD  method  achieves  an  excellent  approximation  of  reflection- 
free  lattice  truncations.  This  is  accomplished  by  the  introduction  of  a 
small,  anisotropic  loss  into  the  region  external  to  the  model  structure  and 
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a wave-field  propagation  condition  at  each  truncation  point.  It  can  be 
shown  that  the  resulting  error  due  to  spurious  wave  reflections  at  the 
lattice  truncations  is  reduced  to  less  than  about  five  percent  by  using 
this  technique. 

3.2.4  The  Plane  Wave  Source  Condition 

Another  basic  problem  with  any  finite-difference  solution  of  Maxwell's 
equations  in  an  unbounded  region  is  the  generation  of  a long  duration  pulse 
or  a continuous  sinusoidal  wave.  Although  such  a wave  could  be  programmed 
as  an  initial  condition,  this  would  result  in  a waste  of  computer  storage 
because  the  lattice  would  have  to  be  elongated  to  initially  contain  the  full 
pulse  or  wave  train.  Another  possibility  would  be  to  vary  the  electi^ic 
field  at  all  points  along  one  end  face  of  the  lattice  in  a pulsed  or  sinu- 
soidal manner.  This  lattice  plane  would  then  radiate  the  desired  plane 
wave  toward  the  model  structure.  However,  such  a specification  of  field 
values  at  a lattice  truncation  plane,  without  consideration  of  the  values  of 
any  possible  outgoing  scattered  waves,  would  cause  undesired  reflection  of 
such  waves  and  significant  error,  as  discussed  in  the  previous  topic. 

The  FD-TD  method  achieves  an  excellent  simulation  of  a long  duration 
pulse  or  a continuous  sinusoidal  wave  without  requiring  any  additional  com- 
puter storage  or  causing  any  additional  wave  reflections.  This  is  accom- 
plished by  using  a program  instruction  which  simulates  the  linear  super- 
position of  an  arbitrary  incoming  plane  wave  with  the  ambient  scattered 
fields  at  all  points  on  a single,  transverse  lattice  plane  located  between 
the  model  structure  and  one  of  the  lattice  truncation  planes.  The  desired 
incident  wave  is  generated  at  the  superposition  plane.  But,  most  impor- 
tantly, any  outgoing  scattered  wave  can  propagate  right  through  the  wave 
source  plane  without  reflection  and  reach  the  lattice  truncation  plane 
beyond  to  be  absorbed.  This  condition  simulates  an  arbitrary  plane  wave 
originating  at  infinity,  and  a scattered  wave  returning  to  infinity,  without 
permitting  any  interaction  between  the  two  waves  except  at  the  model  struc- 
ture. It  can  be  shown  that  the  error  resulting  from  this  simulation  is 
negligible  (less  than  one  percent). 
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3.3  Advantages 

The  FD-TD  method  has  the  following  advantages  relative  to  available 
frequency  domain  and  time  domain  techniques  for  electromagnetic  coupling 
and  shielding  problems. 

1.  The  required  computer  storage  of  the  FD-TD  method  increases  only 
as  the  third-power  of  the  ratio  of  overall  structure  size  to 
spatial  resolution,  regardless  of  the  internal  complexity  of 
the  structure.  Other  computer  techniques  which  require  the 
solution  of  simultaneous  equations  usually  have  a sixth-power 
dependence  on  the  size-to-resolution  ratio  for  complex,  inhomo- 
geneous structures.  This  is  a fundamental  dimensional  advantage 
for  the  FD-TD  method  which  allows  it  to  model  geometries  not 
solvable  by  any  other  procedure. 

2.  The  FD-TD  method  can  model  the  surfaces,  apertures,  and  interiors 
of  complex  structures  in  a straightforward  manner  on  a finite- 
difference  lattice.  Only  a data  card  deck  need  be  punched  to 
specify  the  geometry  of  an  arbitrary  scatterer  or  shield. 

The  maximum  resolution  is  limited  only  by  the  size  of  the 
basic  lattice  unit  cell. 

3.  The  FD-TD  method  can  model  structures  with  square  corners  in  a 
natural  way  if  a cubic  finite-difference  lattice  is  used.  This 
avoids  the  problem  of  dealing  with  the  current  and  field  singu- 
larities often  found  at  corners,  which  lead  to  slow  convergence 
or  erroneous  results  with  other  approaches. 

4.  The  FD-TD  method  allows  a natural  and  simple  treatment  of  the 
following  cases  which  are  difficult  or  impossible  to  handle  in 
any  other  way: 

a.  Dielectrics,  conductors,  and  permeable  media  which  are 
anisotropic  and/or  nonlinear; 

b.  Charged  particles  or  other  ionized  media  within  the 
structure,  due  possibly  to  system  generated  EMP; 
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c.  Irradiation  sources  which  are  either  continuous  wave,  pulsed, 
amplitude  modulated,  or  frequency  modulated; 

d.  Irradiation  source  or  material  boundaries  which  are  moving 
at  relativistic  speeds. 

5.  The  FD-TD  method  allows  a unified  treatment  of  free  space  irradia- 
tion problems,  waveguide  obstacle  problems,  or  combinations  of 
the  above,  i.e.  , waveguide  applicators  for  diathermy  or  hyper- 
thermia, waveguide  antennas,  etc. 

3 . 4 Computational  Details  for  a Uniform,  Cubic  Lattice 

3.4.1  System  of  Finite-Difference  Equations 

Using  the  MKS  system  of  units,  and  assuming  that  the  dielectric  param- 
eters, y,  e,  and  a,  are  independent  of  time,  the  following  system  of  scalar 
equations  is  equivalent  to  Maxwell's  equations  in  the  rectangular  coordinate 
system  (x,  y,  z) : 


- i 

dt  y ' 3z  3y  ^ 


3t 


1 

i ( — L 

y ' 3x  3z  ^ 


an. 


z _ 


3t 


i 

y ^ 3y  3x  ' 


3E 


X _ 


3t 


1 


3E..  , 3H  3H 

= - ( — - - - OE  ) 

3t  e ^ 3z  3x  y 


3t  E ^ 3x  3y  z' 


(la) 

(lb) 

(lc) 

(ld) 

(le) 

(lf) 


Yee^^  originally  introduced  a set  of  finite-difference  equations  for  the 
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system  of  Equations  (la)  - (If).  Following  Yee's  notation,  we  denote  a 
space  point  in  a cubic  lattice  as 

(i.j.k)  = (i6,j6,k6)  (2) 

and  any  function  of  space  and  time  as 

F'^(i,j,k)  = F(i6,j6,k6,n6t) , (3) 

where  6 = 6x  = 6y  = 6z  is  the  space  increment,  6t  is  the  time  increment,  and 
i,  j,  k,  and  n are  integers.  Yee  used  centered  finite-difference  expressions 
for  the  space  and  time  derivatives  that  are  both  simply  programmed  and 
second-order  accurate  in  6 and  in  6t,  respectively; 


3F^i,j.k)  _ F"(i+%,j,k)  - F”(i-%,j,k)  ^ 

3x  6 

n+%  n-% 

aF”(i,j,k)  ^ F(i,j,k)  - F(i,j,k)  ^ Q(g^2j 
3t  6t 

To  achieve  the  accuracy  of  Equation  (4),  and  to  realize  all  of  the  space 
derivatives  of  Equations  (la)  - (If),  Yee  positioned  the  components  of  E 
and  ff  about  a unit  cell  of  the  lattice  as  shown  in  Figure  1.  To  achieve  the 
accuracy  of  Equation  (5),  he  evaluated  F and  H at  alternate  half  time  steps. 
The  result  of  these  assumptions  is  the  following  general  system  of  finite- 
difference  equations  for  the  system  of  Equations  (la)  - (If): 


(4) 

(5) 


n+% 

H^(i  ,j+%,k+%)  = 


6t 

y(i .j+%,k+%)6 


n+% 

Hy(i+%,j ,k+%)  = 


6t 

u(i+%.j ,k+%)6 


n-% 

H^(i  ,j+%,k+%)  + 

""n  n " 

Ey(i ,j+%,k+l ) - Ey(i,j+%,k)  + 


n n 

E^(i,j,k+%)  - E^(i ,j+l ,k+%) 

" n-% 

Hy(i+%,j,k+%)  + 

~n  n 

12(1+1 ,j ,k+%)  - E2(i,j,k+%) 

IE  (i+%,j,k)  - E (i+%,j,k+l) 

1 X X 


1 


(6a) 


(6b) 
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n+% 


6t 

y('i+%.j+%.k)6 


n+1 

Ej^(i+%.j.k)  = 


6t  . 

e(i+%,j ,k)6 


n+1 

Ey(i,j+%,k)  = 


6t 

e(i ,j+%,k)6 


n+1 

E^(i.j.k+%)  = 


<St  . 
c(i ,j,k+%)6 


n-% 

(■>■*■%.  j+%.k)  + 

' n n 

Ej^(i+%J+1  ,k)  - E^(i+%,j,k)  + 

n n 

Ey(i.j+%,k)  - Ey(i+l.j+%.k) 


(6c) 


ri  _ a(i+%»j>k)6t-| 
*-  2e(i+%J.k) -* 

2 e(i  + %Vj  , k)  J 


[1  + 


1 


a(i+%.j  .klTt 
2e(i+%,j.k) 


■] 


n 

Ej^(i+%.j,k)  + 


(6d) 

n+%  n+% 

H2(i+%,j+%.k)  - H^(i+%,j-%,k)  + 

n+%  n+% 

Hy(i+%,j ,k-%)  - Hy(i+%,j ,k+%) 


[1  - 
[1  + 


p(i  ,j+%.k)6t-| 
2e(i.j+%.kT^ 

p(i.j+%.k)6^ 

2e(i  ,j+%;k)  J 


[1  + 


1 

o(i  J+%,k)Tt^ 
2e(i  .j+%,k)  ■' 


n 

Ey(i.j+%.k)  + 


(6e) 

n+%  n+% 

Hj^(i  .j+%.k+%)  - H^(i  ,j+%,k-%)  + 

n+%  n+% 

H2(i-%.j+%.k)  - H^(i+%,j+%,k) 


[1  - 
[1  + 


p(Vj.k+%)5t 

j »k+%) . 
p(i »j .k+%)6t 
2e(i ,j ,k+%) 


I 

-] 


1 

ri  + P(i.j.k+%)6tn 

'■  2e(i  ,j  ,k+%) 


n 

E^(i,j,k+%)  + 


{6f) 

n+%  n+% 

Hy(i+%,j ,k+%)  - Hy(i-%,j ,k+%)  + 

n+%  n+% 

Hj^(i  .j-%.k+%)  - H^(i  ,j+%,k+%) 
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With  the  system  of  Equations  (6a)  - (6f),  the  new  value  of  a field 
vector  component  at  any  lattice  point  depends  only  on  its  previous  value 
and  on  the  previous  values  of  the  components  of  the  other  field  vector  at 
adjacent  points.  Therefore,  at  any  given  time  step,  the  computation  of  a 
field  vector  may  proceed  one  point  at  a time. 

Many  electromagnetic  interaction  problems  involve  nonpermeable  media 
and  can  be  approached  using  a fixed  time  step  and  space  increment.  For  such 
problems  (including  the  cylinder  and  nose  cone  geometries  specified  for  this 
research  effort),  the  quantity  6t/p(i,j,k)6  is  constant  for  all  (i,j,k)  of 
the  lattice,  and  the  Yee  system  of  Equations  (6a)  - (6f)  can  be  simplified 
to  reduce  computer  running  time  in  the  following  manner.  We  define  the 
constants: 

(7a) 
(7b) 
(7c) 

(7d) 

(7e) 

where  m is  an  integer  denoting  a particular  dielectric  or  conducting  medium 
in  the  space  to  be  modeled.  We  also  define  the  proportional  electric-field 
vector 

I = (8) 

Using  the  definitions  of  Equations  (7a)  - (7e)  and  (8),  we  rewrite  Equations 
(6a)  - (6c)  as: 

n+% 

Hj^(i  ,j+%,k+%) 
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n-%  ..n  ^n 

= H^(i ,j+%,k+%)  + Ey(i ,j+%,k+l)  - Ey(i,j+%,k)  + 


^n  ~n 

E (i,j.k+%)  - E^(i ,j+l ,k+%) 


n+% 

n-% 

^n 

^n 

Hy(i+%,j,k+%) 

= Hy(i+%,j,k+%) 

+ E^(i+1 J ,k+%)  - 

E^(i,j,l<+%)  + 

..n 

n 

(9b) 

Ex(i+%.j.k)  - 

E^(i+%.j,k+l) 

n+% 

n-% 

.^n 

= H^(i+%,j+%,k) 

+ Ex(i+%,j+l ,k)  - 

E (i+%,j,k)  + 

~n 

..n 

(9c) 

Ey(iJ+%.k)  - 

Ey(i+1 ,j+%,k) 

This  modification 

eliminates  the  three  multiplications 

needed  by  Yee  in  the 

H part,  of  the  algorithm.  Further,  we  rewrite  Equations  (6d)  - (6f)  as: 
m = MEDIA(i+%,j,k) 


Ej^(i+%,j,k) 


~n 

Cg(m)Ex(i+%,j,k) 


n+%  n+% 

C|j(ni)[H^(i+%,j+%,k)  - H^(i+%,j-%,k)  + 


n+%  n+% 

Hy(i+%.j  .k-%)  - Hy{i+%,j  ,k+%)] 


(9d) 


m = MEDIA{i,j+%,k) 


.^n+1  ^n  n+%  n+% 

= Cg(m)Ey(i  ,j+%,k)  + C|^{m)[Hx(i  ,j+%,k+%)  - Hx(i  ,j+%,k-%)  + 

n+%  n+%  (9e) 

H^{i-%.j+%.k)  - H^(i+%.j+%,k)] 


m = MEDIA{i,j,k+%) 


^n+1  ^n  n+%  n+% 

E^(i.j.k+%)  = Cg(m)E^(i  ,j  ,k+%)  + C|^(m)[Hy  (i+%,j  ,k+%)  - Hy(i-%,j  ,k+%)  + 


n+%  n+% 

Hx(i .j-%»k+%)  - Hx{i ,j+%,k+%)] 


(9f) 


This  modification  eliminates  the  need  for  computer  storage  of  separate  e and 
0 arrays.  Now,  only  a MEDIA  array  which  specifies  the  type-integer  of  the 


dielectric  or  conducting  medium  at  the  location  of  each  electric  field 
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component  in  the  lattice  need  be  stored.  In  addition,  the  and  a of  each 
medium  can  now  be  changed  without  having  to  re-punch  a large  data  card  deck, 
if  the  basic  structure  geometry  is  unchanged.  Such  a change  involves  only 
the  recalculation  of  the  few  values  of  C^(ni)  and  C|^(m). 

3.4.2  Choice  of  Space  and  Time  Increments 

The  choice  of  6 and  6t  is  motivated  by  the  reasons  of  accuracy  and 
algorithm  stability,  respectively.  To  insure  the  accuracy  of  the  computed 
spatial  derivatives  of  the  electromagnetic  fields,  6 must  be  small  compared 
to  a wavelength  (usually  < A/10).  Further,  to  insure  that  the  cubic  lattice 
approximation  to  the  surfaces  of  the  structure  modeled  is  not  too  coarse,  6 
must  be  small  compared  to  the  overall  dimensions  of  the  structure. 

To  insure  the  stability  of  the  time-stepping  algorithm  of  Equations 
(9a)  - (9f),  6t  is  chosen  to  satisfy  the  inequality 


6t  < 


-1 

‘'max 


< 


c /T 
max 


(for  a cubic  lattice) 


(10) 


where  c is  the  maximum  wave  phase  velocity  within  the  model.  The  corres- 
max 

ponding  stability  criterion  set  forth  by  Yee  in  Equations  (7)  and  (8)  of  liis 
paper  is  incorrect.  The  derivation  of  Equation  (10)  is  outlined  as  follows. 

Derivation  of  the  Stability  Criterion 

For  convenience,  we  consider  a normalized  region  of  space  with  p = 1, 
e = 1,0=0,  and  c = 1.  Letting  j = /-T , we  rewrite  Maxwell's  equations  as 

.i  V X (H  + jE)  = 9(H  + jE)/9t,  (lid) 

or  more  simply  as 

j V X V = 'dV/at,  where  V = H + jE.  (lib) 

The  stability  of  a particular  numerical  representation  of  Equation  (11b)  can 
be  examined  simply  by  considering  the  following  pair  of  eigenvalue  problems; 
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(12a) 

^ Inumerical  ^ ^ ~ ' 

(12b) 

Using  the  numerical  time  derivative  given  by  Equation  (5),  Equation  (12a) 
gives 

_n+%  _n-%  _n 

(V  - V )/6t  = XV  . (13) 

_n+%  _n-% 

Defining  a solution  growth  factor  q = V /V  , and  substituting  into 
Equation  (13),  we  solve  for  q: 

q = X6t/2  + /l  + (X6t/2)2.  (14) 

Algorithm  stability  requires  that  |q|  <1  for  all  possible  spatial  modes  in 
the  lattice.  For  this  to  occur. 

Re  X = 0;  | Im  X|  < 2/5t.  (15) 

We  now  let 

V(Jl,m,n)  = V exp[j(k  ll6x  + k m6y  + k n6z)]  (16) 

0 X y z 

represent  an  arbitrary  lattice  spatial  mode.  Using  the  numerical  space 
derivation  formulation  of  Equation  (4),  Equation  (12b)  yields 

sin(%k  6x)  sin(%k  6y)  sin(%k  (Sz)  _ _ 

-2[ ^ ^ ] X V(«,,m,n)  = XV(«,,m,n).  (17) 

After  performing  the  cross  product  and  writing  the  x,  y,  and  z component 

2 

equations,  the  resulting  system  is  solved  for  X : 

o sin^(%k  5x)  sin^(%k  6y)  sin^(%k  6z) 

= -4[ + / + (18) 

6x^  dy*^  6z 

For  all  possible  k , k , k , 

X y z 

Re  X = 0;  | Im  X | < 2(-^- + -^  + -^)  ^ (19) 

6x  6y  6z 

To  satisfy  the  stability  condition  of  Equation  (15)  for  the  arbitrary  lattice 
spatial  mode,  we  set 
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+ 


(20) 


1 


1 % 

— y)  < 2 /6t . 
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The  algorithm  stability  condition  follows  immediately  from  Equation  (20). 

In  an  inhomogeneous  region  of  space,  it  is  difficult  to  determine  a spectrum 
of  X,  analogous  to  Equation  (19),  for  all  possible  lattice  spatial  modes. 

For  absolute  algorithm  stability.  Equation  (10)  suffices  because  it  repre- 
sents a "worst  case"  choice  of  6t. 


3.4.3  Lattice  Truncation  Conditions 

A basic  consideration  with  the  FD-TD  lattice  is  the  treatment  of  the 
field  vector  components  at  the  lattice  truncation  planes.  Inspection  of 
Equations  (9a)  - (9f)  indicates  that  the  values  of  such  components  cannot  be 
determined  from  the  system  of  finite-difference  equations  because  of  the 
centered  nature  of  the  spatial  derivatives.  Therefore,  these  values  must  be 
computed  using  an  auxiliary  truncation  condition.  However,  great  care  must 
be  taken  because  this  condition  must  not  cause  the  spurious  reflection  of 
waves  scattered  outward  from  the  structure  modeled,  as  observed  by  Yee.  The 
goal  of  formulating  the  truncation  condition  is  to  make  the  lattice  trunca- 
tion planes  invisible  to  all  possible  waves  propagating  within  the  lattice, 
as  shown  in  Figure  2. 

A desirable  truncation  condition  relates  in  a simple  way  the  values  of 
the  field  components  at  the  truncation  planes  to  field  component  values  at 
points  one  or  more  6 within  the  lattice.  We  now  consider  examples  of  such 
a truncation  for  cases  of  FD-TD  lattices  in  one  and  three  dimensions. 


One-Dimensiona 1 Case 

For  simplicity,  we  consider  waves  having  only  the  E^  and  components 
and  propagating  in  the  +y  directions.  The  one-dimensional  FD-TD  lattice  is 
simply  a y-directed  line  of  points  having  the  E and  H components  inter- 

Z A 

leaved  and  separated  from  each  other  by  0.5  6y.  The  lattice  is  assumed  to 
extend  from  an  E^  component  at  point  y = 0 to  another  E^  component  at  point 
y = J6y.  A time  step  of  6t  = 6y/c  is  used;  a value  which  is  the  maximum 
allowed  by  the  stability  condition  of  Equation  (10)  for  this  lattice  (6x  = 
6z  = “>). 
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Lattice  truncation  plane 


Fig.  2 IDEAL  FD-TD  LATTICE  TRUNCATION  CONDITIONS 
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Subject  to  these  assumptions,  the  truncation  condition  at  point  y = 0, 
n n-1 

(0)  = (1)  (21a) 

simulates  the  free  space  propagation  of  the  magnitude  of  from  the  point 
•"1"  to  the  truncation  point  "0"  in  one  time  step  (the  free-space  propagation 
delay  implied  by  the  time-step  relation).  This  is  an  exact  truncation  for 
this  lattice  in  that  all  possible  -y-directed  waves  are  absorbed  at  0 with- 
out reflection.  If  we  wish  to  simulate  the  truncation  of  the  lattice  at 
point  y = .16,  the  truncation  condition 

n n-1 

(J)  = E^  (J-1)  (21b) 

is  exact  for  all  possible  +y-directed  waves  at  this  point. 

Three-Dimensional  Case 

Here,  we  consider  waves  having  all  six  field  components  and  propagating 
in  all  possible  directions.  The  lattice  is  assumed  to  extend  from: 

and  components  at  x = %6  to 

and  components  at  x = (I  + %)6; 

E^  and  E^  components  at  y = 0 to 

E^  and  E^  components  at  y = J6; 

E and  E components  at  z = 0 to 

X y 

E and  E components  at  z = K6. 

X y 

A time  step  of  6t  = 6/2c  is  used,  a value  which  is  about  13%  lower  than  the 
maximum  allowed  (6t  = 6//J  c)  by  the  stability  condition  of  Equation  (10)  for 
this  lattice  (6x  = 6y  = 6z  = 6). 

No  simple,  exact  truncation  condition,  analogous  to  Equations  (21a)  and 
(21b),  is  apparent  for  this  three-dimensional  space  lattice.  This  is  because 
we  cannot  assume  the  outgoing  waves  to  be  plane  and  normally  incident  on  one 
lattice  boundary.  At  any  truncation  point,  the  local  angle  of  incidence  of 
these  waves  relative  to  the  truncation  plane  is  unknown.  Further,  several 
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different  waves  having  different  local  angles  of  incidence  may  arrive  at  the 
same  time.  No  simple  truncation  condition  can  account  for  all  of  these 
possibilities.  Therefore,  we  can  arrive  at  only  an  approximate  condition 
that  reduces  the  effective  lattice  boundary  reflection  coefficient  to  an 
acceptable  level . 

A set  of  simple,  approximate  truncation  conditions  that  can  be  used 
with  good  results  is  as  follows. 

At  X = % 6: 
n 

n 


n-2  n-2  n-2  (22a) 

= [Hy  (3/2,j,k-%)  + Hy  (3/2.j,k+%)  + Hy  (3/2 .j  .k+3/2)]/3 

n-2  n-2  n-2 

= [H^  (3/2,j+%,k-l)  + (3/2.j+%,k)  + (3/2.j+%,k+l )]/3 


At  X = (I  + %)  6: 
n 

Hy(I+%,j ,k+%)  = 

n 

H2(I+%,j+%,k)  = 


n-2 


n-2 


n-2 


(22c) 


[Hy  (I-%,j.k-%)  + Hy  (I-%,j,k+%)  + Hy  {I-%.j.k+3/2)]/3 


n-2 


n-2 


n-2 


(22d) 


[H^  (I-%,j+%.k-1)  + (I-%.j+%,k)  + (I-%.j+%,k+l)]/3 


0: 

f^(in,0,k)  = 

n-2 

(i+%,l,k) 

n 

t^(i,0,k+%)  = 

-n-2 

E^  (i,l,k+%) 

At  y = J 6: 

n 

n-2 

E^  (inJ-l,k) 

r^(i,..),k+%)  - 

. n-2 

E^  (i,J-l,k+^) 

(23a) 

(23b) 


(23c) 

(23d) 


PA 


(24a) 


At  z = 0: 

Ex(i+%.j.O)  = 

.n-2 

E"(i.j+%,0)  = 

.n-2 

[Ey  ( 

At  z = K 6: 

n 

Ex(i+%.j.K)  = 

n-2 

Ey(i.j+%.K)  = 

:n-2 

[Ey  ( 

- 1-2 


n-2 

Ex 

n-2 


n-2 

Ex  W 

,n-2 

^ Ey  (i 


n-2 

1 (i+: 

X 

,n-2 


(24b) 

I 

(24c) 

l)]/3 

(24d) 


Equations  (22)  - (24)  allow  the  field  value  at  any  truncation  point  to 
rise  to  approach  the  field  value  of  any  outgoing  wave,  thus  lowering  the 
effective  truncation  plane  reflection  coefficient.  This  is  done  by  modeling 
the  propagation  of  an  outgoing  wave  from  the  lattice  plane  adjacent  to  the 
truncation,  to  the  lattice  plane  at  the  truncation,  in  two  time  steps  (the 
free-space  propagation  delay  implied  by  the  time-step  relation).  The  aver- 
aging process  is  used  to  take  into  account  the  possible  local  angles  of 
incidence  of  the  outgoing  wave  at  the  truncation  and  possible  multiple 
incidences. 

Truncation  conditions  (22)  - (24)  are  useful  for  an  assumed  +y-directed 
incident  plane  wave  with  field  components  and  For  such  a wave.  Equa- 

tion (23)  represents  exact  truncations  similar  to  Equation  (21).  In  addi- 
tion, Equations  (22)  and  (24)  have  no  effect  on  the  propagation  of  such  a 

wave,  which  lacks  H , H , E , and  E . Thus,  this  set  of  truncation  condi- 
y z X y 

tions  effectively  makes  the  lattice  boundary  planes  invisible  to  a +y- 
directed  incident  plane  wave. 

Use  of  Exterior- Region  Anisotropy 

One  way  of  reducing  spurious  reflections  at  the  lattice  truncations  is 
to  introduce  an  anisotropic  lossy  medium  outside  of  the  modeled  structure. 
Properly  constituted,  the  medium  would  attenuate  field  components  present 
only  in  the  scattered  wave,  leaving  the  incident  plane  wave  unaffected. 
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For  the  three-dimensional  case,  this  can  be  easily  done  by  specifyinq  an 

anisotropic  conductivity,  in  the  free-space  region  exterior  to  the 

structure.  Equation  (9f) , the  finite-difference  equation  for  E^,  requires 

= 0 to  insure  that  the  incident  wave  is  not  attenuated.  However,  we 

may  assume  a small  value  of  for  Equations  (9d)  and  (9c),  the  finite- 

difference  equations  for  E and  E , without  affecting  the  propagation  of  the 

A y 

incident  wave  or  the  penetrating  wave  within  the  structure.  This  assumption 

results  in  attenuation  of  the  E and  E components  of  the  exterior  diffracted 

X y 

wave,  and  thus,  reduces  z-directed  wave  reflections  at  the  lattice  trunca- 
tions. 


Effect  on  Algorithm  Stability 

The  stability  condition  of  Equation  (10)  is  valid  for  the  Yee,  or  null 
choice,  of  lattice  truncation  conditions.  This  is  because  Yee's  set  of  trun- 
cation conditions  causes  total  reflection  of  all  lattice  wave  modes  at  the 
surface  planes  of  the  lattice,  and  thus,  introduces  no  new  wave  modes. 
However,  introduction  of  Equations  (22)  - (24)  to  the  three-dimensional 
algorithm  is  found  to  increase  the  strictness  of  the  stability  condition. 

For  the  three-dimensional  case,  some  case  must  be  taken  to  avoid  algorithm 
instability. 

The  nature  of  the  instability  of  the  three-dimensional  algorithm  is  of 
importance.  First,  it  is  late  in  appearance,  requiring  more  than  five- 

4 

hundred  time  steps  for  a 2 x 10  - cell  lattice,  and  more  than  one-thousand 

4 

time  steps  for  a 6 x 10  - cell  lattice.  Second,  its  initial  visibility  is 
delayed  by  either  increasing  the  size  of  the  lattice,  or  by  increasing  the 
losses  of  the  dielectric  media  of  the  lattice.  This  suggests  the  importance 
of  wave  propagation  effects  in  the  growth  of  the  instability. 

There  are  two  likely  solutions  to  the  problem  of  algorithm  instability. 
First,  6t  can  be  reduced.  This,  however,  would  complicate  the  programming 
of  the  truncation  conditions  because  a wave  would  no  longer  propagate  across 
a free-space  unit  cell  of  the  lattice  in  an  integral  number  of  time  steps. 

The  second  solution  is  much  simpler  since  it  does  not  require  pro- 
gramming for  interpolating  field  values  at  the  truncation  planes  between 

5 

time  steps.  This  solution  is  merely  to  set  a lower  bound  of  about  10  cells 
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on  the  size  of  the  FD-TD  lattice  used,  thus  delaying  cne  onset  of  the  insta- 
bility for  several  thousand  time  steps.  In  most  cases,  such  a delay  would 
be  sufficient  for  the  computed  solution  to  reach  the  sinusoidal  steady  state. 
For  many  problems,  the  use  of  a 10  - cell  lattice  is  not  at  all  extravagant 
and  allows  the  problem  of  algorithm  instability  to  be  essentially  forgotten. 
The  computer  runs  in  the  present  research  program  for  the  cylinder  and  the 
missile  nose  cone  employ  this  solution.  They  use,  respectively,  94,000 
cells--800  time  steps  and  58,000  cells--900  time  steps,  without  any  apparent 
instability  of  the  computed  solution. 

3.4.4  Plane  Wave  Source  Condition 

/Mother  basic  consideration  with  the  FD-TD  method  is  the  simulation  of 
the  continuous,  sinusoidal,  incident  plane  wave.  Yee  specified  the  shape  and 
direction  of  propagation  of  an  incident  wave  pulse  by  inserting  all  of  its 
field  values  as  initial  conditions  over  a portion  of  the  lattice.  However, 
the  Yee  approach  is  clearly  inadequate  for  a continuous  wave  train  because  a 
very  elongated  lattice  would  be  needed  to  contain  the  wave  as  an  initial 
condition,  wasting  much  computer  storage. 

In  this  section,  we  discuss  the  simulation  of  an  incident,  +y-directed 
plane  wave  using  a source  condition  localized  at  only  one  lattice  plane,  and 
invisible  to  all  scattered  waves  propagating  within  the  lattice.  This 
allows  a compact  lattice  and  maximum  utilization  of  the  available  computer 
storage. 

The  most  simple  approach  to  this  problem  is  to  vary  the  electric  field 
at  all  points  along  lattice  plane  y = 0 in  a sinusoidal  manner.  This  plane 
would  then  radiate  the  desired  plane  wave.  However,  such  a specification  of 
field  values  at  a lattice  truncation  plane,  without  consideration  of  the 
values  of  the  fields  of  any  possible  outgoing,  scattered  waves,  would  cause 
undesired  wave  reflections. 

A more  desirable  plane  wave  source  condition  would  take  into  account 
the  scattered  fields  at  the  source  plane.  For  the  three-dimensional  case,  a 
useful  wave  source  condition  at  plane  y = 3^6  (near  y = 0)  is  as  follows: 

n ...n 

^ % sin(2Trfn6t)  + E^(i,j5,k+%)  (25) 
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where  f is  the  irradiation  frequency  and  is  defined  by  Equation  (7c). 

Equation  (25)  is  a modification  of  the  Maxwell's  equations  algorithm  for  all 
points  on  the  lattice  plane  y = At  each  point  on  this  source  plane, 

the  computer  first  calculates  E^  in  the  normal  manner  of  the  algorithm,  and 
stores  the  value  in  memory.  Then,  the  value  of  the  sinusoid  is  calculated 
and  added  to  the  stored  value  of  E^.  Finally,  this  modified  value  of  E^  is 
stored  in  memory.  In  effect.  Equation  (25)  simulates  the  linear  superposi- 
tion of  a +y-directed  plane  wave  and  the  ambient  field  along  the  source 
plane.  This  condition  permits  any  scattered,  outgoing  wave  to  propagate 
right  through  the  wave  source  plane  without  reflection,  and  reach  the  lat- 
tice truncation  at  y = 0 to  be  absorbed.  ^ 

3.4.5  Symmetry  Conditions 

An  important  savings  of  computer  memory  and  program  execution  time  ; 

results  if  even  symmetry  of  the  modeled  structure  about  one  or  two  lattice 
planes  can  be  assumed.  In  this  section,  we  discuss  the  programming  of  this 
symmetry  for  the  three-dimensional  case. 

For  the  three-dimensional  case,  the  modeled  structure  is  assumed  to  be 
evenly  symmetric  about  lattice  planes  x = (I  + %)6  and  z = K6; 

e,a( I+%+h,j ,k)  = e,a(I+%-h,j ,k) , P = Pq  (26a) 

e,a(i,j,K+h)  = e,a(i,j,K-h)  (26b) 

The  incident  radiation  is  assumed  to  be  a +y-directed  plane  wave,  with  the 
field  components  E^  and  naturally  having  even  symmetry  about  any  lattice 
plane  x = constant  or  z = constant.  Therefore,  we  conclude  that  the  E^  and 
components  of  the  total  field  possess  even  symmetry  about  the  lattice 
planes  x = (I+%)6  and  z = K6: 


^n  _ ^n 


E2(I+%+h,j ,k+%)  = 

E^( I+%-h,j ,k+%) 

(27a) 

n+% 

(I+%+h,j+%,k+%) 

n+% 

= (I+%-h,j+%,k+%) 

(27b) 

n n 

E2(i.j.K+h)  = E^(i 

J.K-h) 

(27c) 

n+% 

(i,j+%,l<+h)  = 

n+% 

(iiJ+%.K-h) 

(27d) 

28 

To  develop  a convenient  set  of  symmetry  conditions,  we  follow  a procedure 
detailed  below.  This  results  in 

n+l_  n+%_  n+%_ 

(I+%,j+%,l<)  = 0 for  all  n;  (28a) 

n+%  _ _ 

H (i+%,j+%,K)  = E (i+%,j,K)  = E (i,j+%,K)  = 0 for  all  n.  (28b) 

^ ^ y 

Equations  (28a)  and  (28b)  are  sufficient  to  truncate  the  FD-TD  lattice  at 
planes  x = (I+%)6  and  z = K6,  respectively,  by  permitting  the  calculation  of 
the  complete  set  of  field  components,  with  full  specification  of  the  assumed 
even  symmetry. 

Derivation  of  the  Symmetry  Conditions 
A.  Plane  x = (I+%)6 

The  symmetry  conditions  for  the  E , H , and  H field  components  at 
X y z 

lattice  plane  x = (I+%)6  are  derived  from  time-stepping  algorithm  Equations 
(9a)  - (9d)  and  symmetry  assumptions  (26a),  (27a),  and  (27b).  The  first 
step  of  the  derivation  involves  the  determination  of  the  type  of  symmetry 
exhibited  by  E^  about  plane  x = (I+%)6.  To  begin,  we  write  Equation  (9a) 
for  i = I and  for  i = I+l : 


n+%_ 

(I.j+%,k+%) 


n-%_  ^n  _ n _ 

(I.j+%.k+%)  + Ey(I.j+%,k+l)  - Ey(I.j+%,k)  + 


n _ _ 

E^(I,j.k+%)  - E^(I,j+l,k+%) 


(29a) 


n+%_  n-%_  n _ n _ 

(I+l.j+%,k+%)  = H (I+l.j+%,k+%)  + !(I+l.j+%,k+l)  - E„(I+l,j+%,k)  + 
A A y y 


^n  _ ^n  _ 

E^(I+l,j,k+%)  - E^(I+l,j+l,k+%) 


(29b) 


Using  the  symmetry  conditions  of  (27a)  and  (27b)  for  the  case  h = %,  we 
subtract  Equation  (29b)  from  (29a)  and  simplify: 

yT,j+%,k+l)  - E"(T+l.j+%,k+l)  = yT,j+%,k)  - yT+l,j+%,k)  (30) 

Equation  (30)  is  in  the  form  f(j+%,k+l)  = f(j+%,k).  This  implies  a solution 
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is  a constant.  We  thus  have 


of  the  type  f(j+%,k) 


= where 


^n  _ ^n  _ 

Ey(I,j+%,k)  - Ey( I+l ,j+%,k)  = Cj^^,  (31a) 

We  argue  that  the  symmetry  of  the  problem  precludes  the  possibility  of  a con- 
stant step  discontinuity  of  across  the  plane  of  symmetry.  Therefore, 
must  equal  zero  and 

n _ n _ 

Ey(I,j+%,k)  = Ey(I+l,j+%,k).  (31b) 

Equation  (31b)  is  a statement  of  the  even  symmetry  of  E^  about  plane  x = 
(I+%)6. 

Now,  we  may  derive  the  symmetry  conditions  for  E , H , and  H at  plane 
X = (I+%)6.  Using  symmetry  condition  (27a)  for  the  case  h = %,  we  write 
Equation  (9b)  for  i = T: 

n+%_  n-%_  n _ n _ 

Hy  (I+%.j,k+%)  = (I+%,j,k+%)  + E^(I+%,j,k)  - E^(H-%,j,k+l)  (32a) 

Using  the  derived  symmetry  condition  of  Equation  (31b),  we  write  Equation 
(9c)  for  i = T: 

h"  ^+%,j+%,k)  = h"  (T+%,j+%,k)  + E"(T+%,j+l,k)  - E"(T+%,j,k)  (32b) 

From  Equations  (32a)  and  (32b),  we  see  that  Hy  and  at  the  symmetry  plane 
can  depart  from  their  zero  initial  conditions  only  if  E^(T+%,j,k)  assumes 
some  non-zero  values.  Yet,  from  Equation  (9d),  E^  at  the  symmetry  plane  is 
seen  to  remain  at  zero  if  Hy(T+%,j ,k+%)  and  H^(T+%,j+%,k)  are  zero.  There- 
fore, using  an  inductive  argument,  we  conclude  that  these  three  field  com- 
ponents must  remain  at  zero  for  all  time  steps.  Equation  (28a)  is  a state- 
ment of  this  behavior. 

B.  Plane  z = K6 

The  symmetry  conditions  for  the  E , E , and  H field  components  at 
X y z 

lattice  plane  z = K6  are  derived  from  time-stepping  algorithm  Equations  (9c) 

- (9f)  and  symmetry  assumptions  (26b),  (27c),  and  (27d).  The  first  step  of 
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the  derivation  involves  the  determination  of  the  type  of  symmetry  exhibited 
by  about  plane  z = K6.  To  begin,  we  write  Equation  (9f)  for  k = K-1  and 
for  k = K: 

= MEDIA(i  ,j  ,!<-%) 

~n+l  _ _ n+%  _ n+%  _ 

= C^(m^)E^(i  ,j  ,K-%)  + )[Hy(i+%,j  ,K-%)  - ,K-%)  + 

n+%  _ n+%  _ 

H^(i,j-%.K-%)  - H^(i,j+%,K-%)] 


m2  = MEDIA(i ,j ,K+%) 


.n+1 

E^(i.j,K+%) 


n+% 


n+% 


+ Cjj(m2)['^y(i+%«J’K+%)  - H (i-%,j  ,K+%)  + 


n+%  _ n+%  __ 

H^(iJ-%.K+%)  - H^(i,j+%,K+%)] 


(33b) 


Using  the  symmetry  conditions  of  Equations  (26b),  (27c),  and  (27d)  for  the 
case  h = %,  we  subtract  Equation  (33b)  from  (33a)  and  simplify: 

n+%  _ n+%  _ n+%  _ n+%  _ 

Hy(i+%,j,K-%)  - Hy(i+%,j,K+%)  = Hy(i-%.j,K-%)  - Hy(i-%,j,K+%)  (34) 


Equation  (34)  is  in  the  form  f(i+%,j)  = f(i-%,j).  This  implies  a solution 
of  the  type  f(i+%,j)  = C j , where  is  a constant.  We  thus  have 

n+%  _ n+%  _ 

Hy  (i+%,j,K-%)  - Hy  (i+%,j,K+%)  = Cy  (35a) 

We  argue  that  the  symmetry  of  the  problem  precludes  the  possibility  of  a con- 
stant step  discontinuity  of  H across  the  plane  of  symmetry.  Therefore,  C. 

y 3 

must  equal  zero  and 

n+%  _ n+%  _ 

Hy  (i+%,j.K-%)  = Hy  (i+%,j,K+%).  (35b) 

Equation  (35b)  is  a statement  of  the  even  symmetry  of  Hy  about  plane  z = K6. 
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Now,  we  may  derive  the  symmetry  conditions  for  , E , and  H at  plane 
z = U.  Using  the  derived  symmetry  condition  of  Equation  (35b),  we  write 
Equation  {9d)  for  k = K: 

m = MEDIA(i+%.j,K) 

-n+1  _ -'1  _ n+% 

= C^(m)E^(i+%,j,K)  + Cjj(m)[H^(i+%,j+%,K)  - 

n+%  (36a) 

Using  symmetry  condition  (27d)  for  the  case  h - %,  we  write  Equation  (9e)  for 
k = K: 


m = MEDIA(i,j+%,K) 

^n+1  _ _ 

Ey(i,j+%,K)  = Cg(m)Ey(i,j+%,K)  + Cjj(m)[H^(i-%,j+%,K)  - 

n+% 

H2(i+%,j+%,K)] 


(36b) 


From  Equations  (36a)  and  (36b),  we  see  that  E and  E at  the  symmetry  plane 

^ y _ 

can  depart  from  their  zero  initial  conditions  only  if  H^(i+%,j+%,K)  assumes 

some  non-zero  values.  Yet,  from  Equation  (9c),  H at  the  symmetry  plane  is 

seen  to  remain  at  zero  if  E (i+%,j,K)  and  E.(i,j+%,1<)  are  zero.  Therefore, 

X y 

using  an  inductive  argument,  we  conclude  that  these  three  field  components 
must  remain  at  zero  for  all  time  steps.  Equation  (28b)  is  a statement  of 
this  behavior. 


3.5  Review  of  Past  Usages  and  Results 

In  this  section,  we  review  past  usages  and  results  of  the  FD-TD  method 
for  simple  two  and  three-dimensional  dielectric  geometries,  with  an  objective 
of  establishing  the  level  of  accuracy  of  the  method.  Three  sources  of  error 
have  been  considered:  the  approximation  of  the  space  and  time  derivatives  of 
Maxwell's  equations  by  finite-difference  expressions;  the  residual  wave 
reflections  at  the  lattice  truncations;  and  the  stepped-surface  approximation 
of  the  shape  of  a curved  scatterer.  The  computed  results  used  in  this  dis- 
cussion were  obtained  by  the  principal  investigator  during  his  Ph.D.  disser- 
tation work  in  1974  and  1975. 
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3.5.1  Irradiation  of  a Plane  Dielectric  Slab 


First  considered  as  a source  of  error  was  the  approximation  of  the  space 
and  time  derivatives  of  Maxwell's  equations  by  finite-difference  expressions. 
To  isolate  this  error  source,  a scattering  problem  had  to  be  formulated  that 
eliminated  error  due  to  the  approximations  of  the  shape  of  the  scatterer  or 
of  the  lattice  truncations.  Such  a problem  was  the  steady  plane-wave  irradi- 
ation (at  normal  incidence)  of  a flat  dielectric  slab,  because  the  boundaries 
of  the  slab  could  be  defined  exactly  by  two  parallel  lattice  planes.  Further, 
for  a +y-directed  incident  wave,  all  scattered  waves  had  to  propagate  in 
either  the  +y  or  -y  directions;  a situation  where  the  lattice  truncation  con- 
ditions were  exact.  Therefore,  any  observed  error  in  the  results  could  be 
attributed  to  the  finite-difference  approximations  of  the  derivatives. 

For  ease  in  understanding  the  results,  the  slab  was  assumed  to  be  loss- 
less and  one-half  wavelength  thick.  These  conditions  were  fulfilled  by  a 3 
cm  thick  slab,  with  relative  permittivity  equal  to  4,  irradiated  at  2.5  GHz. 
The  geometry  of  this  slab  relative  to  the  problem  lattice  is  detailed  in 
Figure  3a  for  two  lattice  resolutions:  6 = 0.3  cm  = X^/20;  and  5 = 0.6  cm  = 

yio- 

n 

Figure  3b  graphs  the  maximum  absolute  value,  or  envelope,  of  E^{j)  com- 
puted during  the  wave  half-cycle  preceding  the  termination  of  the  algorithm 
n^^  . n , was  chosen  large  enough  to  allow  all  values  of  E within  the 

maX  niaX  Z 

lattice  to  reach  the  sinusoidal  steady  state.  For  the  coarse-lattice  case, 

6t  = 10  psec  and  the  envelope  was  observed  for  130  < n < 150.  For  the  fine- 
lattice  case,  6t  = 5 psec  and  the  envelope  was  observed  for  250  5 n < 300. 

Each  case  represented  an  algorithm  time  of  1500  psec,  or  3.75  wave  cycles  at 
2.5  GHz,  required  to  reach  the  steady  state. 

For  the  fine-lattice  case,  the  envelope  was  virtually  indistinguishable 
from  the  exact  solution  calculated  from  basic  microwave  theory.  In  front  of 
the  slab,  a standing  wave  ratio  of  about  1.04  was  observed.  Within  the  slab, 
the  envelope  decreased  to  exactly  one-half  the  incident  field  magnitude.  In 
back  of  the  slab,  the  envelope  was  constant  and  had  the  same  value  as  the 
incident  field  magnitude.  These  results  were  in  excellent  agreement  with  the 
exact  solution  obtained  either  by  viewing  the  slab  as  a 1:1  impedance 
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Fig.  3 IRRADIATION  OF  A PLANE  DIELECTRIC  SLAB 


transformer  between  its  two  faces,  or  by  using  an  infinite  series  approach. 
The  maximum  error  was  about  +2%.  For  the  coarse-lattice  case,  error  was 
evident  in  the  computation  of  the  envelope  in  front  of  the  slab.  Here,  a 
standing  wave  ratio  of  about  1.13  was  observed.  The  maximum  error  was  esti- 
mated to  be  about  +7%. 

To  insure  that  the  uncertainty  caused  by  the  finite-difference  approxi- 
mations of  the  field  derivatives  is  significantly  less  than  +10%,  it  was  seen 
that  a lattice  resolution  of  6 < A/20  must  be  maintained.  For  inhomogeneous 
scatterers,  A of  this  criterion  should  be  taken  as  the  minimum  value  expected 
within  the  lattice. 

3.5.2  Irradiation  of  a Square  Dielectric  Cylinder 

Next  considered  as  a source  of  error  was  the  presence  of  residual  wave 
reflections  due  to  imperfect  lattice  truncations.  To  isolate  this  error 
source,  a scattering  problem  had  to  be  formulated  that  generated  a roughly 
radially-propagating  diffracted  wave  (to  test  the  truncation  conditions  at 
non-normal  incidence)  without  having  additional  error  due  to  any  stepped-edge 
approximation  of  a curved  boundary  of  the  scatterer.  Such  a problem  was  the 
TM  irradiation  of  an  infinitely- long , rectangular,  dielectric  cylinder,  with 
the  incident  wave  propagating  normally  to  one  cylinder  face.  Here,  the 
boundary  of  the  cylinder  could  be  defined  exactly  by  intersecting  lattice 
planes;  yet,  the  desired  diffracted  wave  would  be  generated.  Any  error  in 
the  results  in  excess  of  that  observed  for  the  plane  dielectric  slab  problem 
could  thus  be  directly  attributed  to  spurious  wave  reflections  at  the  lattice 
truncation  planes. 

For  ease  in  understanding  the  results,  the  cylinder  scattering  problem 

21 

was  assumed  to  have  the  geometry  analyzed  by  Tong  with  a surface  integral 
equation  approach.  Tong's  cylinder  was  composed  of  lossless  dielectric  with 
= 3.84,  and  had  a square  cross  section  with  diameter  d = A^y-ir.  This  con- 
dition was  fulfilled  by  a 3.82  cm  diameter  square  cylinder,  irradiated  at 
2.5  GHz.  The  geometry  of  this  scatterer  relative  to  the  problem  lattice  is 
detailed  in  Figure  4a  for  two  lattice  resolutions:  6 = 0.095  cm  = A^/64, 
for  cylinder  perimeter  ABCD;  and  6 = 0.191  cm  = A^/32,  for  cylinder  perimeter 
A'B'C'D' . 
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Figure  4b  graphs  the  computed  results  for  the  envelope  of  the  cylinder's 
surface  electric  field,  E^(s),  along  perimeters  ABCD  and  A'B'C'D',  as  well 
as  Tong's  solution.  (Here,  s represents  a normalized  position  along  the 
perimeter.)  For  the  ABCD  solution,  6t  = 1.59  psec  and  the  envelope  was 
observed  for  630  5 n < 756.  For  the  A'B'C'D'  solution,  6t  = 3.18  psec  and 
the  envelope  was  observed  for  315  < n < 378.  Each  case  represented  an  algor- 
ithm time  of  1200  psec,  or  3.0  wave  cycles  at  2.5  GHz,  required  to  reach  the 
steady  state. 

Since  the  lattice  resolutions  for  both  cases  of  Figure  4 were  finer  than 
that  of  the  slab  problem  of  Figure  3,  it  was  inferred  that  the  error  due  to 
the  finite-difference  approximations  of  the  derivatives  was  less  than  +2% 
for  each  case.  Any  error  above  this  limit  was  assumed  due  to  residual  wave 
reflections  at  the  lattice  truncations.  Over  most  of  the  surface  of  the 
cylinder,  the  error  level  of  each  FD-TD  solution  was  comparable  and  limited 
to  about  +10%.  The  principal  disagreement  between  the  FD-TD  solutions,  and 
the  largest  error  of  the  ABCD  case,  occurred  at  the  electric  field  minimum. 
Evidently,  the  interaction  of  the  cylinder  and  the  lattice  truncations  was 
so  weak  that  it  could  influence  the  computed  results  only  a field  minima. 

It  was  concluded  that  the  FD-TD  lattice  truncation  conditions  lead  to  an 
error  level  of  less  than  +10%  at  most  points,  even  for  small  spacings 
between  the  structure  modeled  and  the  lattice  truncations. 

3.5.3  Irradiation  of  a Circular  Dielectric  Cylinder 

The  stepped-surface , or  lattice,  approximation  of  the  shape  of  a curved 
scatterer  was  next  considered  as  a source  of  error.  To  estimate  this  error 
source,  the  TM  irradiation  of  an  infinitely-long,  circular  dielectric 
cylinder  was  modeled.  Here,  in  addition  to  error  due  to  imperfect  lattice 
truncations,  modeling  error  was  introduced  because  the  cylinder  surface 
could  not  be  exactly  defined  by  combinations  of  intersecting  lattice  planes. 
In  effect,  spurious  wave  diffraction  effects  were  caused  at  each  intersec- 
tion of  two  planes.  The  magnitude  of  this  additional  error  could  be  esti- 
mated by  comparing  the  accuracy  of  the  circular  cylinder  model  to  that  of 
the  square  cylinder  model  of  Figure  4. 
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The  geometry  of  the  circular  cylinder  model  relative  to  the  problem 
lattice  is  detailed  in  Figure  5a.  The  lattice  coordinates  internal  to  the 
cylinder,  determined  by  [(i-24.5)  + j-25)  ]'^  < 20,  were  assigned  the  dielec- 

tric parameters  = 4.0  and  a = 0.  An  irradiation  frequency  of  2.5  GHz  was 
assumed,  with  6 = 0.3  cm  = X^/20,  6t  = 5 psec,  and  the  envelope  of 
observed  for  460  5 n < 500  (a  maximum  algorithm  time  of  2500  psec,  or  6.25 
wave  cycles).  The  cylinder  thus  had  a diameter  of  1 free-space  wavelength. 

Figures  5b  and  5c  detail  the  computed  values  of  |E  (24,j)|/E7.  and 

z ^inc 

|E7(15,j)|/E2.  , respectively.  These  figures  also  present  the  exact  solu- 

tion  calculated  using  the  summed-mode  series  technique  of  Jones.  As  seen 
in  the  figures,  the  FD-TD  solution  located  the  positions  of  the  peaks  and 
nulls  of  the  electric  field  with  a maximum  error  of  +6,  or  about  +3%  of  the 
diameter  of  the  cylinder.  The  magnitude  of  each  peak  was  determined  with  a 
maximum  error  of  +10%.  Error  in  excess  of  10%  appeared  in  the  computation 
of  the  magnitude  of  several  of  the  field  minima.  Overall,  the  uncertainty 
of  this  solution  was  found  to  be  comparable  to  that  of  the  square  cylinder 
solution.  Evidently,  the  error  caused  by  the  stepped-surface  approximation 
of  the  circular  cylinder  was  small  in  comparison  with  that  caused  by  the 
lattice  truncation  conditions. 

3.5.4  Irradiation  of  a Dielectric  Sphere 

To  further  investigate  the  error  introduced  by  the  lattice  approximation 
of  a curved  scatterer,  the  plane-wave  irradiation  of  a dielectric  sphere  was 
modeled.  Unlike  the  circular  cylinder,  the  sphere  has  a surface  which  must 
be  approximated  in  three  dimensions.  The  magnitude  of  any  additional  error 
caused  by  this  approximation  could  be  estimated  by  a simple  comparison  of 
the  sphere  results  with  those  of  the  circular  cylinder. 

The  geometry  of  the  sphere  model  relative  to  the  problem  lattice  is 
depicted  in  Figure  6 at  the  two  lattice  symmetry  planes.  The  lattice  coor- 
dinates  internal  to  the  sphere  were  determined  by  [(i-19.5)  + (j-20)  + 

(k-19)  < 15.  To  allow  a close  comparison  with  the  circular  cylinder 

results,  all  of  the  dielectric  and  program  parameters  of  that  run  were 
repeated,  namely;  = 4.0;  a = 0;  f = 2.5  GHz;  6 = 0.3  cm  = Xj/20;  6t  = 5 
psec;  and  envelope  observation  for  460  < n < 500.  The  sphere  thus  had  a 
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diameter  of  0.75  free-space  wavelength.  To  reduce  spurious  reflections  at 
the  lattice  truncations  (as  discussed  in  Section  3.4.3),  a value  of 
equal  to  0.1  mho/m  was  assumed. 

Figures  7a  and  7b  detail  the  computed,  normalized  values  of  two  electric 

field  components  near  the  sphere  irradiation  axis:  | E (19. 5 ,j  ,18) | /E^ • 

j 1 n c 

and  |E  (19,j,18.5)|/E7.  . These  figures  also  present  the  exact  solution 

z ^ 1 n c o o 

calculated  using  the  summed-mode  series  technique  of  Stratton.^  As  seen  in 
the  figures,  the  FD-TD  solution  located  the  positions  of  the  peaks  and  nulls 
of  the  electric  field  with  a maximum  error  of  +6,  or  about  +2%  of  the  diam- 
eter of  the  sphere.  The  magnitude  of  each  peak  was  determined  with  a maximum 
error  of  +10%.  Overall,  the  uncertainty  of  this  solution  was  found  to  be 
comparable  to  that  of  the  square  and  circular  cylinder  solutions.  Evidently, 
the  error  caused  by  the  stepped  approximation  of  the  surface  of  the  sphere 
in  three  dimensions  was  small  in  comparison  with  that  caused  by  the  lattice 
truncation  conditions.  It  was,  therefore,  concluded  that  the  three  dimen- 
sional FD-TD  program  allows  solutions  with  a level  of  accuracy  comparable  to 
that  of  the  two  dimensional  program. 
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4.0  DESCRIPTION  OF  THE  PRESENT  RESEARCH  PROGRAM 


4 . 1 Introduction 

Electromagnetic  coupling  and  shielding  problems  have  traditionally  been 
difficult  to  treat  with  analytical  or  numerical  methods  because  of  the  fail- 
ure of  these  methods  to  adequately  resolve  the  effects  of  the  apertures, 
curvatures,  corners,  and  internal  contents  of  structures.  A practical 
numerical  method  has  been  desired  to  allow  the  direct  modeling  and  solution 
of  complex  realistic  coupling  problems.  The  purpose  of  this  research  program 
was  to  evaluate  the  suitability  of  the  finite-difference  time-domain  (FD-TD) 
solution  method  for  Maxwell's  equations  to  determine  the  amount  of  electro- 
magnetic coupling  through  an  aperture  into  an  enclosed  conducting  container. 

The  FD-TD  method  allows,  in  principle,  the  computation  of  the  internal 
fields  of  complex  conducting  geometries.  However,  before  the  present 
research  program,  this  method  had  not  been  utilized  and  evaluated  for  any 
conducting  geometries.  To  build  up  confidence  in  the  FD-TD  method  for 
future  applications,  it  was  desired  to  evaluate  the  usage  of  the  method  for 
certain  simple,  generic  metal  structures. 

During  the  present  research  program,  two  specific  metal  structures  were 
used  in  this  evaluation:  the  first,  an  aluminum  cylinder  with  one  open  end; 
the  second,  the  nose  cone  section  of  a missile.  Each  structure  was  modeled 
using  the  FD-TD  method  to  compute  the  internal  EM  fields  generated  by  an 
incident  plane  wave  propagating  along  the  structure  axis.  The  results  of 
the  cylinder  model  were  then  compared  to  available  theoretical  and  experi- 
mental data.  Final  evaluation  of  the  nose  cone  results  will  be  possible 
when  reliable  experimental  data  for  this  geometry  is  obtained  in  a future 
research  program. 

The  following  section  describes  each  coupling-analysis  problem  con- 
sidered and  summarizes  the  steps  taken  to  solve  each  problem.  Detailed 
descriptions  of  the  results  follow  in  succeeding  sections. 


4.2  Description  of  C 


s Problems  Considered 


4.2.1  Task  1:  Prediction  of  the  Coupling  Into 
an  Open-Ended  Cylinder 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  - Circular  (19.0  cm  diameter),  68.5  cm  long, 
open-ended  aluminum  cylinder,  as  shown  in  Figure  8; 

Incident  wave  - 300  MHz  plane  wave  propagating  down  the  cylinder 
axis  toward  its  open  end; 

Desired  fields  - Each  component  of  total  E and  total  H in  the  axial 
cross-section  plane  of  the  cylinder  down  to  40  cm  from  the  open  end. 
First,  with  the  cross-section  plane  parallel  to  the  incident  E,  and 
again  with  the  plane  parallel  to  the  incident  H; 

Resolution  - 0.5  cm  uniformly  throughout  the  mapping  planes; 

Plotted  values  - In  decibels  relative  to  an  incident  E of  1 volt/ 
meter  and  an  incident  H of  1/377  ampere/meter. 


To  solve  this  coupling  problem,  an  existing  FD-TD  computer  code  was 
suitably  modified  to  model  the  cylinder  of  Figure  8.  The  following  steps 
were  taken: 


a.  The  existing  FD-TD  computer  code  was  modified  for  the  Control 
Data  STAR-100  computer.  A 24  x 163  x 24  - cell  lattice  was 
programmed,  with  even  symmetry  of  the  incident  fields  and 
cylinder  assumed  about  lattice  planes  x = 24.56  and  z = 246. 

b.  The  cylinder  geometry  of  Figure  8 was  mapped  into  the  new 
finite-difference  lattice  for  a unit  cell  diameter  of  6 = 0.5  cm. 

c.  The  FD-TD  program  was  run  for  this  geometry  for  800  time  steps 
(equivalent  to  2.0  cycles  of  the  incident  wave)  assuming  loss- 
less air  within  the  cylinder.  10^  words  of  memory  and  3.5 
minutes  of  central  processor  time  were  required  on  the  STAR-100. 
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OPEN-ENDED  ALUMINUM  CYLINDER  GEOMETRY  OF  TASK  1 (not  to  scale) 


d.  The  computed  results  for  E and  H along  the  cylinder  axis  were 
plotted  at  intervals  of  200  time  steps  (0.5  cycle  of  the  inci- 
dent wave).  Slow  convergence  of  the  results  to  the  sinusoidal 
steady  state  was  observed. 

e.  The  FD-TD  program  was  re-run,  this  time  assuming  a slight 

amount  of  isotropic  loss  = 0.01  mho/m)  for  the  air  with- 

in the  cylinder  to  increase  the  rate  of  convergence  to  the 
steady  state.  Again,  800  time  steps  were  completed. 

f.  Plotting  of  E and  H along  the  cylinder  axis  for  the  second  run 
indicated  a much  more  rapid  convergence  to  the  steady  state. 

The  final  computed  results  for  the  fields  were  then  reduced  to 
contour  maps  along  the  symmetry  planes,  and  to  simple  graphs 
at  selected  cuts  through  the  symmetry  planes. 

4.2.2  Task  2:  Prediction  of  the  Coupling  Into 
a Missile  Nose  Cone 

The  FD-TD  technique  was  employed  to  solve  the  following  electromagnetic 
coupling  problem: 

Interacting  structure  - Aluminum  nose  cone,  shown  in  Figure  9,  with 
two  apertures:  a circular  one  in  the  nose,  and  a sleeve  fitting 
located  23-1/3  cm  aft.  Missile  body  geometry  beyond  sleeve  fitting 
assumed  to  continue  to  infinity  with  constant  cross-section  shape. 
Aperture  cases  investigated: 

Trial  1 - Sleeve  fitting  open,  nose  aperture  closed. 

Trial  2 - Both  apertures  open; 

Incident  wave  - 300  MFIz  plane  wave  propagating  down  the  axis  of  the 
structure  toward  its  nose  aperture; 

Desired  fields  - Each  component  of  total  E and  total  FI  in  the  axial 
cross-section  plane  of  the  nose  cone.  First,  with  the  cross-section 
plane  parallel  to  the  incident  E,  and  again  with  the  plane  parallel 
to  the  incident  F; 
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Fig.  9 NOSE  CONE  GEOMETRY  OF  TASK  2 (not  to  scale) 


Resolution  - 1/3  cm  uniformly  throughout  the  mapping  planes; 


Plotted  values  - In  decibels  relative  to  an  incident  E of  1 volt/ 

meter  and  an  incident  H of  1/377  ampere/meter. 

To  solve  this  coupling  problem,  the  following  steps  were  taken: 

a.  The  24  x 163  x 24  - cell  lattice  of  the  aluminum  cylinder  pro- 
gram was  truncated  to  24  x 100  x 24  cells. 

b.  The  geometry  of  Figure  9 was  mapped  into  the  new  finite-difference 

lattice  for  a unit  cell  diameter  of  6 = 1/3  cm. 

c.  The  FD-TD  program  was  run  for  the  aperture  case  of  Trial  1 for 

900  time  steps  (equivalent  to  1.5  cycles  of  the  incident  wave), 
assuming  a slight  amount  of  isotropic  loss  = 0.025  mho/m) 

within  the  cylinder  to  speed  the  rate  of  convergence.  6 x 10^ 
words  of  memory  and  2.8  minutes  of  central  processor  time  were 
required. 

d.  The  FD-TD  program  was  run  for  the  aperture  case  of  Trial  2 for 
900  time  steps,  assuming  the  same  interior  loss, 

e.  For  each  trial,  the  final  computed  results  for  E and  H were 
reduced  to  contour  maps  along  the  symmetry  planes. 

4.3  Details  and  Results  of  Task  1: 

Coupling  Into  an  Open-Ended  Cylinder 

4.3.1  Cylinder  Model 

Figures  10  and  11  depict  the  geometry  of  the  cylinder  model  used  for 
Task  1.  Figure  10  shows  the  stepped-surface  model  of  the  cylinder  wall  used 
for  lattice  planes  j = 14  through  j = 150.  This  model  was  specified  by  set- 
ting the  conductivity,  a,  equal  to  that  of  aluminum  (3.7  x 10^  mho/m)  for 

individual  E , E , and  E components  nearest  the  desired  circular  locus. 

X y z 

This  resulted  in  the  modeling  of  the  cylinder  wall  as  an  aluminum  sheet 
having  virtually  zero  thickness.  Figure  11  shows  the  positioning  of  the 
cylinder  relative  to  the  front  and  back  planes  of  the  lattice.  As  shown  in 
this  figure,  the  cylinder  aperture  was  located  at  plane  j = 14,  and  the 
cylinder  back-plane  was  located  at  plane  j = 151.  To  the  front  and  rear  of 
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Fig.  10  CROSS-SECTION  OF  CYLINDER  MODEL  FOR  TASK  1 
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Fig.  11  VIEW  OF  CYLINDER  MODEL  FOR  TASK  1 
AT  HORIZONTAL  SYMMETRY  PLANE 
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the  cylinder  lay  regions  of  air.  In  the  front  air  region,  the  simulated  300 
MHz  incident  wave  was  generated  at  plane  j -■  3 with  the  field  components 
and  H . 

X 

For  all  lattice  cells  exterior  to  the  model  cylinder,  the  anisotropic 
conductivity,  equal  to  0.01  mho/m  was  assumed  to  help  improve  the 

lattice  truncation  conditions,  as  discussed  in  Section  3.4.3.  This  value  of 
Oext  caused  the  exponential  decay  of  and  E^  fields  in  the  exterior  region 
with  an  effective  skin  depth  of  about  110  lattice  cells. 

Using  the  computer  program  listed  in  Appendix  A,  the  cylinder  model  was 
completely  specified  by  punching  3 groups  of  24  cards,  giving  a total  of  72 
cards.  Card  group  1 specified  the  air  medium  of  the  lattice  in  front  and 
back  of  the  cylinder  (planes  J = 0 through  j = 13  and  planes  J = 152  through 
j = 163).  Card  group  2 specified  the  stepped-surface  model  of  the  cylinder 
(planes  j = 14  through  j = 150).  Finally,  card  group  3 specified  the 
cylinder  backplane  (plane  j = 151).  The  format  of  the  data  cards  is  dis- 
cussed in  Appendix  A. 

4.3.2  Convergence  of  the  Computer  Fields 

Two  800  time-step  programs  (each  equivalent  to  2.0  cycles  of  the  inci- 
dent wave)  were  run  during  Task  1 in  order  to  investigate  the  rate  of  con- 
vergence of  the  computed  fields  to  the  sinusoidal  steady  state.  The  first 
program  assumed  lossless  air  within  the  cylinder;  the  second  program  assumed 
a small  isotropic  conductivity,  equal  to  0.01  mho/m  within  the  cylinder. 

The  purpose  of  modeling  the  slightly  lossy  air  was  to  cause  the  reactive 
fields  within  the  cylinder  to  converge  more  rapidly  to  the  expected  beyond- 
cutoff  condition. 

Figures  12  and  13  are  graphs  of  the  computed  total  electric  field, 

|E,/E7.  I (in  decibels)  along  the  axis  of  the  cylinder  for  the  cases  of 
^int  ” ^ °int  ~ mho/m,  respectively.  In  Figure  12,  curves  are 
plotted  for  the  cases  n = 400  and  n = 800  time  steps;  in  Figu  e 13,  curves 
are  plotted  for  n = 200,  n = 400,  n = 600,  and  n = 800  time  steps.  Each 
curve  gives  the  computed  field  envelope  during  the  200  time-steps  period 
(0.5  cycle  of  the  incident  wave)  before  the  specified  value  of  n.  In  Figure 
13,  it  should  be  noted  that,  at  200  time  steps,  the  incident  wave  penetruted 
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Distance  from  aperture  (cm) 


Fig.  12  COMPUTED  ELECTRIC  FIELD 
ALONG  CYLINDER  AXIS 
(a.  ^ = 0 CASE) 


only  about  45  cm  into  the  cylinder,  causing  the  sharp  downward  break  of  the 
n = 200  curve.  For  all  other  curves,  the  incident  wave  penetrated  fully  to 
the  cylinder  backplane. 

Comparing  Figures  12  and  13,  it  is  seen  that  the  use  of  caused  the 
linear  decibel  slope  at  the  cylinder  aperture  to  lengthen  from  a maximum 
depth  of  12.5  cm  to  more  than  25  cm  after  800  time  steps.  Further,  the  ulti- 
mate computed  field  attenuation  increased  frcm  about  30  dB  to  almost  55  dB. 
This  improvement  in  the  convergence  of  the  cylinder's  internal  fields  to  the 
expected  cutoff  condition  was  much  more  than  the  11  dB  of  field  decay  that 
o^.^^  would  cause  for  a wave  propagating  the  full  length  of  the  cylinder. 
Evidently,  the  transient  internal  fields  were  highly  reactive  (carried  little 
real  power  flow)  and  quickly  dissipated  when  forced  to  supply  energy  to  main- 
tain an  electric  field  distribution  across  a slightly  lossy  medium. 

Figure  13  also  indicates  the  rate  of  field  convergence  to  the  steady 
state  under  the  condition  of  finite  Once  the  wave  Tielcis  were  estab- 

lished throughout  the  cylinder  (curves  n = 400,  600,  and  800),  the  principal 
effect  of  an  added  200  program  time  steps  was  the  lengthening  of  the  linear 
decibel  slope  extending  from  the  aperture,  and  a consequent  deepending  of 
the  field  null  within  the  cylinder.  Noting  the  break  points  of  the  curves 
from  the  linear  decibel  slope,  each  200  time-step  increase  is  seen  to  have 
lengthened  the  slope  by  about  10  cm  and  decreased  the  residual  internal 
fields  by  about  10  dB. 

4.3.3  Comparison  With  Available  Data 

Electric  Field  Along  the  Cylinder  Axis 

Using  results  of  the  FD-TD  program.  Figure  14  graphs  the  computed 

IE  /E,.  I along  the  cylinder  axis  as  a solid  curve.  The  FD-TD  results  are 
' z ^inc' 

after  800  time  steps  with  o-  ^ = 0.01  mho/m  (the  n = 800  curve  of  Figure  13). 

I nt  p. 

For  comparison,  experimental  results  are  shown  as  circled  points,  and  com- 

24 

puted  results  using  the  BOR-3  body  of  revolution  code  are  graphed  as  a 
dashed  curve. 

From  Figure  14,  it  is  seen  that  the  FD-TD  curve  lies  between  the  BOR-3 
theory  and  the  experimental  results.  The  FD-TD  and  BOR-3  curves  have  nearly 


54 


the  same  slope  down  to  about  -55  dB,  where  the  FD-TD  curve  levels  off. 

However,  the  experimental  results  level  out  at  only  -30  dB.  A previous 

study  of  the  experimental  procedure  has  shown  that  this  saturation  was  due 

to  inadequate  suppression  of  radiated  power  at  the  third  harmonic  of  the 

nominal  test  frequency,  which  excited  an  above-cutoff  mode  in  the  cylinder.^^ 

This  study  has  also  indicated  that  the  experimental  results  were  consistently 

1-2  dB  above  the  likely  actual  values  due  to  problems  in  the  design  and 

25 

mounting  of  the  test  probe.  Consideration  of  these  experimental  uncer- 
tainties improves  the  correlation  between  the  test  data  and  the  FD-TD  results. 

Overall,  the  FD-TD  method  is  seen  to  yield  results  for  |E  /E7.  I along 

z ^inc 

the  cylinder  axis  which  agree  well  with  available  ,oretical  and  experi- 
mental results  in  terms  of  both  relative  decay  rate  and  absolute  magnitudes. 
Next,  comparison  of  data  will  be  made  for  the  electric  field  along  vertical 
and  horizontal  radial  lines  of  the  cylinder. 

Electric  Field  Along  Vertical  Radial  Lines 

Using  results  of  the  FD-TD  program.  Figure  15  graphs  the  computed 

lE^/E,.  I along  vertical,  radial  lines  of  the  cylinder  (parallel  to  E7.  ). 

z ^inc'  j ^inc 

Separate  solid  curves  are  plotted  for  radial  lines  at  distances,  d,  of  0 cm, 
10  cm,  20  cm,  and  35  cm  from  the  aperture.  Again,  the  FD-TD  results  are 
after  800  time  steps  with  = 0.01  mho/m.  For  comparison,  experimental 
results  are  shown  as  circled  points  for  the  cases  d = 0 cm  and  d = 10  cm. 
(Experimental  results  at  greater  depths  are  not  shown  because  of  their  rapid 
saturation  at  -30  dB  due  to  the  third-harmonic  test  problem.) 

From  Figure  15,  it  is  seen  that  the  experimental  results  at  d = 0 cm 
and  d = 10  cm  are  all  within  2 dB  of  the  FD-TD  results.  This  correlation  is 
very  good;  yet  it  might  be  improved  upon  consideration  of  the  likely  1-2  dB 
"high"  readings  reported  for  the  test  probe. 

Electric  Field  Along  Horizontal  Radial  Lines 

Figure  16  is  similar  to  Figure  15  except  that  I 1 graphed 
along  several  horizontal  radial  lines  of  the  cylinder  (parallel  to 
The  experimental  results  are  seen  to  be  within  3 dB  of  the  FD-TD  results 
over  most  of  the  range  of  the  results. 
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Fig.  15  COMPUTED  ELECTRIC  FIELD  ALONG  VERTICAL 
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Fig.  16  COMPUTED  ELECTRIC  FIELD  ALONG  HORIZONTAL 
RADIAL  LINES  (PARALLEL  TO  ) OF 

THE  CYLINDER 
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4.3.4  Computed  Field  Maps 

Using  results  of  the  FD-TD  program  (800  time  steps,  = 0.01  mho/m). 

Figures  17  and  18  graph  contour  maps  of  the  computed  field  components  at  the 
cylinder's  vertical  and  horizontal  symmetry  planes,  respectively.  All  com- 
ponent magnitudes  are  normalized  to  either  IE7.  I or  iHy.  i = IL,.  1/377, 
and  are  given  as  decibel  numbers.  Contours  are  plotted  at  exact  6 dB  inter- 
vals by  using  a linear  interpolation  method  to  determine  each  contour's 
position  between  adjacent  field  envelope  points.  Although  the  lattice  cell 
diameter,  6,  equals  0.5  cm,  this  interpolation  method  allows  the  generation 
of  smooth  curves  in  most  cases  without  a 0.5  cm-period  stair-case  effect. 

It  should  be  noted  that  only  three  field  components  are  non-zero  at 

each  symmetry  plane:  , and  E^  at  the  vertical  plane;  and  E^,  H^,  and 

H at  the  horizontal  plane.  This  Tact  was  dei ived  in  Section  3.4.5  of  tnis 
y 

report. 
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(a)Ez  (b)Hx  (c)Ey 

Fig.  17  FD-TD  FIELD  CONTOURS  IN  CYLINDER  VERTICAL  SYMMETRY  PLANE 
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Distance  From  Open  End  (Cm) 
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Distance  From  Axis  (Cm)  Distonce  From  Axis  (Cm)  Distonce  From  Axis  (Cm) 

(a)E2  (b)  (c)  Hy 

Fig.  18  FD-TD  FIELD  CONTOURS  IN  CYLINDER  HORIZONTAL  SYMMETRY  PLANE 
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4.4  Retails  and  Results  of  Task  2: 

Coupling  Into  a Missile  Nose  Cone 

4.4.1  Nose  Cone  Model 

Figures  19  and  20  depict  the  geometry  of  the  nose  cone  model  used  for 
Task  2.  Figure  19  shows  the  stepped-surface  model  of  the  nose  cone  wall 
used  at  the  front  circular  aperture  and  at  the  rear  sleeve-fitting  aperture. 

A value  of  a equal  to  3.7  x 10^  mho/m  (aluminum)  was  assigned  to  the  field 
components  at  the  wall.  Figure  20  shows  the  positioning  of  the  nose  cone 
relative  to  the  front  and  back  planes  of  the  lattice.  As  shown  in  this 
figure,  the  nose  aperture  was  located  at  plane  j = 11,  and  the  sleeve  fitting 
was  located  at  planes  j = 74  - ^2.  To  the  front  of  the  nose  cone  lay  a 
region  of  air;  to  the  rear  of  the  nose  cone,  the  missile  body  was  assumed 
to  extend  to  infinity  with  a constant  circular  cross  section.  Based  upon 
the  FD-TD  results  of  the  cylinder  problem  of  Task  1,  this  infinitely-long- 
missile  assumption  was  not  expected  to  cause  significant  error.  (In  fact, 
very  little  wave  reflection  had  been  computed  at  the  cylinder  backplane  at 
intermediate  time  steps.)  In  the  front  air  region,  the  simulated  300  MHz 
incident  wave  was  generated  at  plane  j = 3 with  the  field  components  and 
H . 

X 

For  all  lattice  cells  exterior  to  the  model  nose  cone,  the  anisotropic 
conductivity,  equal  to  0.025  mho/m  was  assumed  to  help  improve  the 

lattice  truncation  conditions,  as  discussed  in  Section  3.4.3.  This  value  of 
a ^ caused  the  exponential  decay  of  E and  E fields  in  the  exterior  region 
with  an  effective  skin  depth  of  about  65  lattice  cells.  To  speed  convergence 
of  the  interior  fields,  as  discussed  in  Section  4.3.2,  the  isotropic  con- 
ductivity, equal  to  0.025  mho/m  was  selected  for  the  cylinder  interior. 

Using  the  computer  program  listed  in  Appendix  A,  the  nose  cone 
geometries  for  both  Trial  1 (sleeve-fitting  open,  nose  aperture  closed)  and 
Trial  2 (both  apertures  open)  were  specified  by  a total  of  464  data  cards. 
Only  16  cards  h^d  to  be  re-punched  to  change  from  Trial  1 to  Trial  2. 
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Fig.  19  KEY  CROSS-SECTIONS  OF  NOSE-CONE  MODEL  FOR  TASK  2 


missile 


Fig.  20  VIEW  OF  NOSE  CONE  MODEL  FOR  TASK  2 AT 
HORIZONTAL  SYMMETRY  PLANE 


4.4.2  Computed  Field  Maps 

Trial  1 --  Only  Sleeve-Fitting  Open 

Using  results  of  the  FD-TD  program  (900  time  steps,  = 0.025  mho/in), 

Figures  21  and  22  graph  contour  maps  of  the  computed  field  components  at 
the  nose  cone  vertical  and  horizontal  symmetry  planes,  respectively,  for 
Trial  1.  Contours  are  plotted  at  exact  10  dB  intervals  using  linear  inter- 
polation. The  intersections  of  the  darker  grid  lines  (spaced  by  5 minor 
divisions,  a total  of  6 = 1/3  cm)  denotes  the  location  of  the  field  vector 
components  in  the  lattice  symmetry  planes.  Because  of  the  staggered  posi- 
tions of  the  field  components  around  a lattice  unit  cell,  these  darker  grid 
lines  may  vary  in  position  relative  to  the  fixed  nose  cone  walls  by  +"/2  = 
+1/6  cm,  or  2.5  minor  divisions. 

It  should  be  noted  that  the  stepped-surface  approximation  of  the  smooth, 
tapered  nose  cone  wall  introduces  cusp-like  distortions  in  seve.'ol  of  the 
field  contours.  However,  these  distortions  are  only  manifested  within  about 
1 cm  of  the  point  of  each  surface  step.  Very  likely,  the  exact  field  contour 
here  can  be  found  simply  by  drawing  a smooth  curve  connecting  the  adjacent 
undisturbed  contour  sections. 

Trial  2 --  Both  Apertures  Open 

Again,  the  FD-TD  program  was  run  for  900  time  steps  with  = 0.025 
mho/m.  Figures  23  and  24  graph  contour  maps  analogous  to  those  of  Figures 
21  and  22.  Comparing  the  corresponding  maps  for  Trial  1 and  Trial  2,  it  is 
seen  that  opening  the  nose  aperture  had  very  little  effect  upon  the  field 
contours  near  the  sleeve  fitting.  Coupling  between  the  two  apertures 
occurred  only  at  field  levels  lower  than  -40  dB.  For  this  reason,  it  was 
decided  not  to  run  a Trial  3 (only  nose  aperture  open)  because  the  resulting 
field  contours  near  the  nose  would  almost  certainly  be  the  same  as  for  Trial 
2. 
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(o)  (b)  (c)  Ey 

Fig.  21  FD-TD  FIELD  CONTOURS  IN  NOSE  CONE  VERTICAL  SYMMETRY 
PIANE.  TRIAL  1-  ONLY  SLEEVE- FITTING  APERTURE  OPEN 
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Distance  From  Axis  (Cm)  Distance  From  Axis  (Cm)  Distance  From  Axis  (Cm) 

(0)62  (b)  Hx  (c)Hy 

Fig.  22  FD-TD  FIELD  CONTOURS  IN  NOSE  CONE  HORIZONTAL  SYMMETRY 
PLANE.  TRIAL  1-  ONLY  SLEEVE -FITTING  APERTURE  OPEN 
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Distonce  From  Axis  (Cm)  Distance  From  Axis  (Cm)  Distance  From  Axis  (Cm) 

(a) Ej  (b)  (c)  Hy 

Fig.  24  FD-TD  FIELD  CONTOURS  IN  NOSE  CONE  HORIZONTAL  SYMMETRY 
PLANE.  TRIAL  2-  BOTH  SLEEVE  AND  NOSE  APERTURES  OPEN 


5.0  DISCUSSION  AND  CONCLUSIONS 


This  research  program  demonstrated  that  the  FD-TD  method  can  be 
successfully  applied  to  axial-incidence,  electromagnetic  coupling  problems 
involving  highly  conducting  structures  with  hole  and  sleeve-type  apertures. 
Accuracy  of  the  FD-TD  results  was  very  good  relative  to  the  uncertainties  of 
available  experimental  and  numerical -theory  approaches.  Convergence  of  the 
EM  fields  to  the  sinusoidal  steady  state  occurred  within  about  2 cycles  of 
the  incident  wave  when  a slight  was  assigned  to  the  structure  interior. 
This  resulted  in  program  running  times  of  3.5  minutes  or  less  on  the 
CDC  STAR-100  for  10^  - cell  lattices. 

The  FD-TD  method  appears  to  have  great  promise  for  applications  involv- 
ing conducting  structures  at  arbitrary  angles  of  incidence,  and  combination 
conducting  - dielectric  structures.  Its  ability  to  achieve  finely-detailed 
models  of  the  interiors  of  such  structures  could  be  utilized  to  determine  the 
internal  fields  of  many  practical  objects.  The  incorporation  of  potential 
computer  program  accelerators,  such  as  variable  cell  size,  could  result  in 
even  more  cost-effective  results. 
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APPENDIX  A 

FD-TD  COMPUTER  PROGRAMS 


A.l  INTRODUCTION 

This  appendix  documents  the  computer  programs  written  during  the  present 
research  effort.  Included  are  listings  of  the  programs  for  the  cylinder 
problem  and  the  nose  cone  problem  and  a description  of  the  data  card  format 
employed.  All  computer  programs  were  written  using  STAR  Fortran  Version  2.1 
for  processing  by  the  Control  Data  STAR-100  computer  system  under  its  1.2 
operating  system.  This  Fortran  version  contains  certain  extensions  to  stan- 
26 

dard  Fortran  that  permit  usage  of  the  vector  processing  capabilities  of 

the  STAR-100.  The  reader  is  referred  to  the  STAR  Fortran  Manual  for 

27 

detailed  discussion  of  these  features. 

A . 2 PROGRAM  LISTING  FOR  TASK  1:  COUPLING  INTO  AN  OPEN-ENDED  CYLINDER 

The  following  9 pages  list  the  computer  program  for  the  24  x 163  x 24 
cell  --  800  time  step  run  of  the  cylinder  problem  (Task  1).  In  the  listed 
problem  parameters,  FREQ  = 3.0  E + 8 denotes  the  operating  frequency,  f = 

300  MHz;  DX  = 0.005  denotes  the  lattice  cell  diameter,  6 = 0.005  m;  MPR 
denotes  the  total  number  of  media  within  the  model,  equal  to  3 (lossless 
air,  aluminum,  slightly  lossy  air);  DATA  EPS  and  DATA  SIG  give  the  assumed 
relative  dielectric  constant  and  conductivity  (mhos/m)  of  each  medium;  and 
MINDT  and  MAXDT  give  the  number  of  the  first  and  last  time  step  of  the 
al gori thm. 


A-1 


o O 


TTPT  68.5  CM  LONG,  OPEN-tN^EO  CONDUCTING  CYLINDER 


TWCrPCNT  WAVE  HflS  THC  COMPONENTS  Er~T^n 'FrX,'~flN[J  IS  DIRECTt: 
ALONG  THE  CYLINDER  AXIS  INTO  ITS  OPEN  END 

7^  X ~163  X 24  CELL  CUtfTC  LATtTCE  rs~ DSED 

UNIT  CELL  diameter  = Ox  = 0.5  CM  = WA VELENGTH/200 


-Ts- 


»OX  IS  ASSUMED 


SOFT  TEM  wave  source  CONDITiaXTAT  PCANE  Y ' = IS  LTSEl 

SOFT  LATTICE  TRUNCATIONS  ARE  USED 

PROGRAM  IS  OPTIKTZED  FOR  ~rKEr-COC  SrA-R-lTTO 


REAL  MUR,  

DIMENSION  A<97350<»)  .Z(5936),Al(60  0),A2{m0  0),AAA(600), 


T aA  ( 25 ) , N(25)'  ,mT(7i,)  ,Or(2A»  ,CAm)  , C B ( 9 ) ,trRSTTrf STTrm 

COMMON  A 

DESCRIPTOR  Dl  ,D7>~D3 ,D» , 05  tNED ,~3V~,  6w , BX  « >^  Y , A 1 X , A2X  ,~A  l Y rATT i 


1 A?r 

BIT  BV,BBB(600) ,BW,B8W (600) , 3X , 8B X ( 600 ) ,BY ,R0 Y ( 2A ) 

ASSIUN  A 1 XTff 1( ZTTSTTT— 

ASSIGN  A2X,A2(27;573) 

ASSIGN  AlY,AH26»^An 

ASSIGN  A2Y»A2(26»574) 

ASSIGN  A1Z*A1 ( 1 ;599)  ~ 

assign  A2Zf A2(1 1599) 

AsS  I GN"  { 1 (TT 

ASSIGN  BX,BBX(i;600> 

assign  By,yBY(II24J  “ 

assign  DA,D0( 1 ;24) 

ASSIGN  D5,DE  (1124) 

T1  = SECOND(CP) 

PRINT  ISO,  T1 

150  FORMAT  {F20.5) 


C I.  PROBLEM  PARAMETERS 

FREU  = J.0Er<^ 

OX  = 0.005 
MPR  ^ J' 

DATA  EPS  / l.O,  l.Ot  1.0  / 

DATA  51G  } Onn  3T7F*7"J  OTOT  7 

MINOT  = I 

HAXDr=  BOTJ 

C 

C II.  BASIC  AND  DLRIy/ED  CONSIANTS 

PI  = 3.14159265 

HUZ  = A. a '*  PI  * l.ftE-7 

EPSZ  = 8.054E-12 
OT  = DX  7 6.0E*B 
MHALF  = 0.5  / FREQ  / OT 


Aid  »600»  = 0. 


C PLAMES  14  - 150  (UnI^’OPM  CIRCULAR  CYLINDER) 

»r«o  4<  f fd  r » t-t  ♦5^ro^  ■ 

WF  AO  <*«  (/ ( I ) « I =^60l « 3200) 

9r«o  T7rTT.rtT^rr»t^?tmT 

t I > I o * I so 

T Tj ' • d-rr  

« # . I'  ^ I T /(  I ISO  is) 

• ^ •IAC'PL  A»*f  


,-•■  » a f 

r'  " 


c 

c 


c 

c 

c 


c 

c 


r 

c 


c 

c 


c 


c 

c 


PRINT  150*  T2 

TV.--  TlHE-STgPPrNTS"  roCTT^ 

00  200  J=MrN0T»MAXDT 

TERW  = SINIFLOAT  (J)  *R0) 

MCAL  = 3 ♦ IFIXIFLOAT ( J)/6.0) 
HCALL  s 

IF (MCAL.LT.54)MCALL=MCAL 


transverse  plane  =l 

rx-i"  EZ  TRUNCATIONST 

A(626I625)  = A(n625) 

■=  A(6S62r5?51 
A(3226;625)  = A(260U625) 


EY  ITERATION...... 

ASSIGN  Olt.OrN.574 

ASSIGN  02»  .OyN.-STA: 

ASSIGN  03*.DyN.574 

A(1926)  = 0.5  *■  < A (2551 ) ♦A7?S52TI 

A(1927I22)  = 0.333  ♦ ( A ( 255 U 22 ) ♦ A ( 2552 ; 22 ) ♦ A ( 2553 » ?2 ) ) 

a 71 9AVT'^~(n  333333  * (A^(25  7J)  ^“?V(r»AT?57^m 

A{2551 t2A)  = A( 1951 ;24) 

Ul“=  CA  TJ)  * A { 195115  K) 

02  = A(3926;574)  - A(39ei;574)  ♦ A(5276;57«»)  - A(5277»57A) 

TJT“"=  UB  TTr-^»T12 

A(1951»574)  = D1  ♦ 03 

ATI'^ZISTETJ  0 1 = QHVCTRCrrOTDTBXt  A { r926  ? CT  D ) ) 

free 


TRANSVERSE  PLANES  ? - 163 

iRj  ^2  ■jri^ri'HCALL — ■ 

JOEL  = (JY-1)  * 17808 


EX  ITERATION... 

TS'S^GN  UI».0YN.5/3 

ASSIGN  02».DYN.573 

ASSIGN  NF0*.0YN.573 
ASSIGN  BV.BBBT 13573 1 

DD"3U  HA=1»3 

M = JOEL  ♦ 5936*(MA-i) 


SOFT  lattice  truncation 

A(65b3*H»  = 0.5  » ( A rTl8H*M)  ♦A  (71B9*Mn 

A(656<**n;21  ) = 0.333  • ( A ( 7 1 88*M » 21 ) ♦ A ( 71  89*M  ; 21 ) ♦ A ( 7 190  ♦M 
I 121 ) ) 

~Tlf65B5*Ml  ^-0.5  ■ (A(7?09<W>  ♦yr7210*Hn 

A(718H*mI23)  ■ AI6588*MI23) 


C 


c 

MAIN  EX  LOOPS 

NFD  = A (59b J) 

DO  ?2  JJ=1*MPR 

k5V  * NhU«t.U*JJ 

AIX  = Q8VCTRL(CA( JJ) fSvlAlX) 

id 

AZX  = UHVLIKL (Cb ( JJ)  »«VIA2X> 

01  = AIX  * A(6588*M|573) 

D2  = Adltfl  J*HI57J)-A  (52  7 1051  J*MI5n)-A(10b38*H 

1 »573) 

03  = A2X  * 02 

■I 

C 

c 

EY  ITERATION 

ASSIGN  Ulf. Lrm.^7ii^ 

ASSIGN  02*.0YN,57A 

ASSIGN  l)3»»UYN»57^ 
assign  NF0».0YN.574 

C 

Assign  bv. bhbtitsta > 

Do  40  MA=1.3 

M s:  JOEL  ♦ 5936*(MA-1) 

C 

C 

SOFT  LATTICE  TRUNCATION.,... 

ATTee^^fTJ  = u.5^  rSlOtfi3^*Mj4.A(845  5*m  ) 

A(7863^m;22)  = 0.333  ® ( A ( 84S7*M 1 22 ) ♦ A ( » 22 ) ♦ A ( 8489*M 
1 »?2)) 

A (T8S5*M)  = 0 • 333  * ( A ( 8509*9 ) • 2t0  * A ( H51 0 ♦*^7') 
A(8487*M;2A)  = A(7887»M;24) 

C 

MAIN  EY  LOOPS 

DO  32  JJ=l»MPR 

BV  *“  NrO#EQ*<JJ 

AlY  = Q8VCTRL (CA ( JJ) ,8V;a1Y) 

AZT  = UBTC:  rTTU  TCm  JJI  .OV.AdTl 

01  = AlY  • A(7887>M»S74) 

- 

Od  - Ai'^8o2*Mt874)“Ai9B3^*M*5/‘*>»Aiil21d’*w«574j”Att3  213*M, 

1 STAY 

03  = A2Y  * 02 

c 

A t ^ OB  7 ♦ H 1 5 1** ) = 01  ♦ Oj 

c 

.....EY  ENVELOPE  COMPtTT  A T I ON  . , , , , 

04  = Q8VCMPRS ( A ( 7861 ♦M»b00) » 3X 104) 

ITS — = v«D^rLnr«  u^i — 

BY  = 05.GT.A (0512*MI^4) 

A (85I2*F*i:4»  = IJHVL  1 HL  (U5*MT  ; A (a512*"»  24)  1 

free 

C~ 

C 

EZ  ITERATION 

assign  OI..OTV.S9N 

ASSIGN  02».0YN.S99 


c 

-ASSIGN  03i.0YN.599 

ASSIGN  NF0».0YN.599 
-ASSIGN  ev-*9B8^t  *59») — 

M = JOEL  ♦ 5936*(MA-1) 

c 

MAIN  EZ  LOOPS 

tHf  u - « 

DO  42  JJ=1,MPR 

BV  - NFO»EQ«JJ 

, AIZ  = 08VCTRL (CA ( JJ) ,BV»A12) 

«2Z  = UMVLTHL (CO fJJ) *HV) A2Z) 

01  = AIZ  * A (9162*M)599) 

1 599) 

03  = A2Z  * 02 

C 

^rr9  T 6 2 ♦ MT5  99)  = TJl  ♦ U3 

L 

r7TTTE7"3Dn  ihM  irfAvt  souRCt  toNu I f roTj;";".  .'i 

IF ( JY.GE.2.0R.MA.LE.2)G0  TO  47 

C 

A (21034)599)  s TERM  ♦ A ( 2 1034 ) 599 ) 

C 

97 

....  .TZ  envelope  COMPUTA TION. . . . . 

04  = VABSIA (9737*M)24) )D4) 

-Err  -=:-D4.GT  . A (9  ^8  T*H;24) 

A(9787*M)24)  = Q0VCTWL(O4.8Y;A(9787*M)24) ) 

U4  = UHVLMKKb  ( A (9  Ibl  » -iX  » U4  » 

05  = VAbS(D4»05) 

50 

BY  = 05. G I .A  (9017^)24) 

A(9012*M)24)  = Q8VCTRL (D5*8Y) A (9012*M)24) ) 

C 

C 

HX  ITERATION 

DO  60  MA=l«3 

M = JOEL  ♦ 5936»(MA-1) 

C 

C 

main  HX  LOOPS 

□1  = Jrrj9U 1 ♦M5599) 

02  = A ( 1951 *4)599) -A< 1926*m;599) ♦A (3226*m;599) -A{9162*M» 

1 599) 

C 

A (390r*M)5g9)  = Dr^*TJ2 

L 

. . . . .HX  ENVELOPE  L’OHPUTAT ION.  , ... 

04  = VABS( A (4476*M)24) )D4) 

BY  = AT4B2b*^rZ4r^ 

A(4526*M;24)  = Q8VCTRL(04»BY)A(4526*M«2A) ) 

04^' = OHVCWFBSnrr3V0(r»BT600)  *3X104) 

05  = VA8S(D4)05) 

60 

BY^  a 13B.UT  .A  r455T^*Tr243 

A(4551*m)24)  « Q«VCTMLi05.8Y)A(4S51*MI24) ) 

FOFI 

C 


AA(2I22)  = 0.333  * ( A ( 1 1 1 37^M I 22) ♦ A ( 1 1 1 38« M 1 22) ♦ A ( 1 1 1 39*M; 
I 2?)  ) 


AAA(i;600)  = Q8VXPNn(AA(i;24) «3WIAAA(1;600) ) 


c 

A(  iii3T*h;24)  = aavcw>^ffS(AriqsT3»mgt)0)  .Bw;ATni3r*'m24) ) 

c 

MAIN  HY  LOOPS..... 

01  = A ( 10513*M»598) 

D2  = A(yibJ*Mlby8i-A(yib2*H»59H)  ♦Ar6b^X*M»'59aJ-Atb^s8»*M» 

1 S9FT 

A{10513»M;598)  = Dl  ♦ D2 

C 

A( 10512*MI600)  = Q8VCTRL(AAA(1;600) ,Bw;a(10S12*M»600) ) 

C 

C 

MT  tNVtLUPt  

04  = VABS(A(11087*HJ24) 104) 

70 

HY  = 04.  t»  I . A ( 1 1 1 6^*M  » 24  J 

A(11162*MI24)  = Q8VCT«HD4»By»A(lll62*Mj24)  ) 

C 

"n?Ei 

C 

.....HZ  1 TLPA T I UN « . . . . 
assign  0l«.0YN.573 

C 

assign  02V.T)VN.57J 

DO  80  MA=If3 

M = JOEL  ♦ 5936* (MA-1)  4 

C 

C 

SOFT  LATTICE  TRUNCATION 

AA(2)  = 0.5  ^ raf5877*M) ♦A(5p79*M>  > 

AA(3«21)  = 0.333  • (A(5877*M;21 ) ♦A(sa78*M«21 ) *A (5879*Mj?i )) 

AA<2"»)  — 0.5  » <A<5o'fa*Hj  • <it5899»M)) 

AAA(H600)  = Q8VXPN0(AA  < H24)  ,3W  JAAA  ( 1 ;t,00)  ) 

C 

A (bo  = UaVUHKK’b  \ 52bZ^  tOW*  A 1 T ♦ dWTl 

C 

HZ  LOOPS#*### 

01  = A(5277*HJ573) 

A (5277*M»573)  =01  * 02 

A (S251  •MI600  > « Q8VCTRL  < AAA  ( I t 600  ) .RWI  A {S?S1  •Mtf.OO)  ) 

“»a  » •»  o»vrr»»t  tiy.t>»^»rAr^2^»*»T^frtrrr 


QO 


free 

-82-corfriNtrr 


IF (MCALL.LE.53)G0  TO  94 


C 

c 

transverse  Plane  ei64 

c 

M = 162  • 5936 

c — 

•••••EXf  EZ  T RUNC A T I ONS • • • • • 

A (6562*M;625)  = A (5937*M;625) 

« 1 bv J f ^ w ^ r dJ  f j 

A (7237 ♦Ml  625)  = A (626* Ml 625) 

A (9162*M»625J  = A ( 0537*M 16253 

A (8537*m;625)  = A ( 9837*M ; 625 ) 

C 

JR 9 8 3T^W»b25I  L32^D ♦ «lb25) 

.....HA  IltRAliuN..... 

ASSIGN  D1*.0YN.599 

ASSiljN  U2*.UTN.6yN 

01  = A(3901*m;599) 

I 599^r 

A (3901 ♦M;599)  = 01  ♦ 02 

04  - yABS<A(4476*M<24)  104) 

BY  = 04.GT.A<4526*M;24) 

A(45g6*H;24)  = QflVCTKL  (D4»BT<  A-t45g6*MI247  7 

04  = Q8VCMPRS( A (3900*M(600) *9X;04) 

05  - VABS(04-r05T 

BY  = D5,GT.A(4551*M;24) 

“A  (4551^1^24)  -=  Q8VCTm.  r05»BYI7n-4 551*m;24)  T 

. free 


C HZ  ITERATION..... 

ASSIGN  D1»»DTN.573 

assign  D2t.DYN.573 

AA(2)  = 0.5  • (A  (5877*Mr»A(5879»»1)  ) 

AA(3J21)  = 0.333  * ( A ( 5877^M ; 21 ) ♦ A ( 5878*M » 2l ) ♦ A ( 5879*M ; 21 ) ) 

AA(241  = 0.5  » (A(bH98*H)  » A(SH99*Hn 

AAAdlSOO)  = Q8VXPND(AA(1I24)  .aw;AAA(l»600)  ) 

A(b876»W;24l  = Q8VCHPR5rA(5252*M;600) «BW;A(S87S»NI24)  ) 

Dl  = A (5277*M|573) 

U2  = A(bbH8*»«IS73>-A(652*H;573)»An95~r^IS73)~A(1952*H;s73) 

A (5277*M;573)  = 01  ♦ 02 

ATS'^bTVRIbUO  J = TJBVCIRLTAAAd  »bOO)  .BW I A ( S25 1 tM  1 500  ) » 

A (5251 ♦MI600)  = Q8VCTRL ( 0 . 0 . BX » A ( 5251 ♦M 1 600 ) ) 

fitee 

94  T3  = SECONO(CP) 

PPim  ISO*  TJ  ~ 


»‘IELD  ENVELOPE  ROUTINE 

00  100  L ■‘‘‘^AlF  .M**OT  .^'HAi  f 
rr  ( J.fo.l  I"*  > TO  101 
100  ' ■»•*<  i Sul 


lo  u 


c 

C AT  horizontal  SYMMETRY  PLANE 

nn~PRTNT“  1 ■ 

102  FORMAT  ( IHl  »52X,27HEZ  ENVELOP*^  FOR  TIME  STEP  =»IS* 

I //.62X,15HPLANE  Z = ?a»D)<  , // ,2X . IH J, /) 

CALL  ENV (3850* 1 .0E*4/RB) 

FR INT  lOJf  (N(LI) »LI  = i,24) 

103  FORMAT (//»8X»24I5»2X» iHl) 


PRINT  10^»  J 

lO'VFORHAI  ( lHi-,!3'2X»2^HH)ri:TIVFLUPF  r UR  TIME  ?IE^^  =»Ib»' 
1 //.62X. 15HPLANE  Z = ?A»DX * // ♦ 2X , 1 H J, / ) 

CALL  ENV(A-525»3.77Eri'67" 

PRINT  103»  (N(LI) ♦LI=1»2A) 


PRINT  105»  J 

— 10b  F0MHAniHl*52X,2/HHY  ENVELOPE  *•  OR  IIME  STEP  =«Ib* 
1 //.62X, 15HPLANE  Z = 2A»0X » // . 2X , IH J , / » 

-CALX  TNV  r5e2S»3.7TE* 61 

PRINT  103.  (N<LI ) .LI=1 .2A) 


AT  VERTICAL  SYMMETRY  PLANE 

PRINI  106.  J 

106  FORMAT ( IHl .52X.27HEZ  ENVELOPE  FOR  TIME  sTE»  =.I5. 

1 /7T62ir,  ITHPLANr  X = ■24.5*0XV77V?X,1TTJ»71 

CALL  ENV <3875. 1 .0£*A/R8) 

PT7TNT  107.  tN(LI)TCT^72AT — 

107  FORMAT (//.8X. 2415. 2X. IHK) 

c 

PRINT  108.  J 

njB  EORflAK  IHl. b2XY2/HHx  ENVELOPE  ^OR  TIME  STEP'  =.'T5. 

1 //.62X. 17HPLANE  X = 24. 5*0X . // . 2X . 1 H J. / ) 

CALL  ENV(45b0.3.7/E*6» 

PRINT  107.  (N(LI ) .LI=1 .24) 

c 

PRINT  109.  J 

— 109  FURMAt (IHl. 52X.27HEY  ENVELOPE  fqr  r i ME  SlEP  =. 15. 

1 //.62X.  17HPLANE  X = 24. 5‘»0X  . // . 2X  . 1 H J.  / ) 

CALL“  ENV  <2575. 1 .0E*4/RB) 

PRINT  107.  <N<LI ) .LI=1 .24) 


199  CONTINUE 

— ZOU  CONTTMOE 

T4  = SECONO<CP) 
PRINI  150.'  F4 
STOP 


— SUqffOirTTNC  evv-rLOC;sCAL^) 

dimension  IPR (25) .a (973504) 

COMMON  A 

00  2 LA=l,162 

LB  s 163  - LA 

LC  = L6  ♦ 1 

HJ-  =-  LB«»S^6-  ♦ -ttJC 

lPR(i;24)  = SCALE  * A(l0»1;24) 

— A7rD*l»24)  - 0. 

2 PRINT  3»  LC.  (IPR(LF) ♦LF=1,24) 

3 FORMAT  ( IX. 13. 5X.24-r5) 

return 

“"END 


A. 3 PROGRAM  LISTING  FOR  TASK  2:  COUPLING  INTO  THE  NOSE 
CONE  OF  A MISSILE 

Most  of  the  Fortran  statements  of  the  nose  cone  program  are  identical 
to  those  of  the  cylinder  program  of  the  preceding  appendix  section.  The 
necessary  changes  involve  only  the  reduction  of  the  length  of  the  main 
data  storage  vector.  A,  and  new  data-read  cards.  Therefore,  in  this  section 
we  list  only  the  modifications  of  the  cylinder  program. 


PROGRAM  FOTO  < I NPUT , OUTPUT  * T APE60= I NPi JT ) 


C 

C PUN  B-  STEADY  300  MhZ  TEH  IRRADIATION  OF  A 12.0  CH 

C oiahttep*  aluminum  nusf  Cone 

C trial  2-  both  sleeve  fitting  and  nose  apertures  are  open 

:: INCIDENT  WAVE  HAS  ThE  COMPONENTS  E7  AND  HX«  AND  IS  DIRECTE 

C ALONG  THE  NOSE  CONE  AXIS 

C 2^  X 100  X 24  CELL  CUBIC  LATTICE  IS  USED 

C UNIT  CELL  diameter  = Ox  = 0.3,1  CM  = WA VELENOTH/lOO 


DIMENSION  A (599536)  , Z ( 5936 ) » A I ( 600 ) , A2 ( 600 ) . A AA ( 600 ) . 

J AA(25).  N(25) *00(24) *0E(2A) .CA(9) .CB(9) .EPS (3) »STG(3) 


C I.  PROBLEM  parameters 

FREQ  = 3.0E*e 


Ox  = 0.01/3.0 

MPR  = 1 

DATA  EPS  / 1.0* 

.1.0* 

1.0  / 

3-7E*7* 

mindt  = 1 

MAXDT  = 9QQ 

III.  LOAD  VECTOR  A. 

ZERO  INITIAL  FIELDS 

Z ( 1 ;5936)  - 0. 

A(i;5936)  = 0. 

A (593601)5936)  = 0. 


TYPE  OF  medium.. 

READ  4,  (Z  ( I ) .1  = 1*600) 

4 FORMAT  (75F1.0) 

READ  4*  (Z(  I )*  1 = 2601  *3200) 
OFflO  4*  (Z ( I ) * 1 = 1 301  * IROO) 
no  5 1*2*10 

inri  ■*  ( i-i ) • 

5 A ( lOF  l • I « Z(|IS9)^) 

».  I A • I I • «>  t . ? 

*1  At  • • / ' I • 1 ' I • 4 (.l  w 

• • i I . , I . 1. 


> C't 


READ  (Z(I) ♦I=I301f 1900) 
lOEL  = (lA-n  ♦ 5936 
A ( I0EL*1 ;5936)  = Z(i;5q36) 
READ  4,  (Zd)  f I = l»600) 

REAP  4»  (Z(I> ♦I=2601<3?00) 

00  6 10=1.6 

lOEL  = (IA»IB-1)  * 5936 

6 A ( IDEL*1 ;5936)  = Z(i;5q36) 

■ DO  7 1=09.100 

lOEL  = (I-l)  * 5936 

7 A ( IOEL*l  <5936)  = Z(H5q36') 

REAP  4.  (Z  ( I ) .1  = 1.600) 

READ  4.  (Z(I) .1=2601. 3?01) 

READ  4.  (7  ( 1 ) .1  = 1301 .IQJ;) 

DO  0 1=74.00 

IDEL  = (I-l)  « 5936 

0 A ( IDEL>1 ;5936)  = Z(i;5936) 

DO  9 1=03.07 

IDEL  = ( I-l ) * 5936 
■9  _A(I0EL*i;5936)  = Z(i;59361 

READ  _4^  _(Z(JJ  .1=1301  . IQQO) 

I = 01 

LIIEJ = (I-l)  « 5936 

A (IDEL* 1 55936)  = Z(i;5936) 


REAP  4.  ( Z ( I ) . I = 1 .600 ) 

RF..An  -4t-  17  111  .T=?60]  t3?QQl 
REAP  4,  (Z ( 1 ).  1 = 1301 , 1900) 

I = 82 

IDEL  = (I-l)  5V36 

A ( IDEL* 1 55936)  = Z(i;5936) 

READ  (Z (I ) . 1 = 1 .600) 

REAP  4,  (Z ( I ). 1=2601 .3200) 
READ  4.  (Z ( I ). 1=1301 . 1900) 
I = 00 

IDEL  = (I-l)  » 5936 

A ( rOEL» 1 55936)  = Z(155936) 


T2  = SECOND(CP) 

— .P.&LMI_15.Q.  T2  

I1M£-5TEPP1MG  LOOP 

00  200  J***1N0T  .»*A»OT 
Tfum  « s IN (FLOAT (J) •MU)  • PP 
^ \ . IT  1 1 (T  I OAT  ( ji 

1»  all  . . ) M Av_>  - 3 » 


ft2  CONTINUE 

IF(MCALL.LE»32)GQ  TO  94  • 

C 

C TRANSVERSE  PLANE  =101 

M = 99  * 5936 


End 


SUBROUTINE  ENV ( LOC , SCAlE ) 

DIMENSION  IPR(g5)«A(59 9536) 

COMMON  A 

DO  2 LA=I,99 

LB  = 100  - LA 

LC  = LB  ♦ 1 

LD  = LB*593'6  ♦ LOC 

IPR  ( 1 ;?4)  = SCAI  F A (lD»I  ;?4) 

A(LD*1I24)  = 0. 

2 PRINT  3,  LCt  (IPR(LF) fLF=I .?4) 

3 F0RMAT(1X,I3,5X»24IS) 

RFTURN 

END 


) 


A. 4 DATA  CARD  FORMAT 


The  data  cards  specify  the  type  of  medium  at  each  location  of  an  elec- 
tric field  component.  Up  to  9 distinct  media  can  be  specified  within  a 
lattice. 


Using  the  Fortran  statement  4 FORMAT  (75F1.0)  a medium-type  integer 
1,  2,...,  9 can  be  assigned  to  the  600  locations  of  an  electric  field  com- 
ponent in  one  plane  j = constant  with  only  8 data  cards.  The  600  locations 
are  ordered  as  shown  in  Figure  A-1.  With  the  75F1.0  format,  we  have 

Data  Card  Assigns  Type  Integer  to 

Consecutive  Locations 


1 

2 

3 

4 

5 

6 

7 

8 


I - 75 

76  - 150 

151  - 225 

226  - 300 

301  - 375 

376  - 450 

451  - 525 

526  - 600 


In  all  data  cards,  column  25,  column  50,  and  columns  75  - 80  are  left  blank 


FIGURE  A-1  ORDERING  OF  THE  ELECTRIC  FIELD  COMPONENT  LOCATIONS 
IN  THE  LATTICE  PLANE  j = CONSTANT 


The  media-typ 
one  plane  j = cons 
with  plane  j = J-i 


56728 


1.8 


ita  for  each  field  component  are  read  into  the  program 
; at  a time,  beginning  with  plane  j = 1 and  ending 
I the  following  order: 


Field 

Component 


Number  of 
Data  Cards 


8 

8 

8 


8 

8 

8 


i 


E 

E 

E 


X 

z 

y 


8 

8 

8 


groups  can  be  used  if  the  system  geometry  is  inde- 
iber  of  lattice  planes.  For  example,  see  Fortran 
DO  5,  DO  6,  and  DO  7 of  the  cylinder  program  (Run  A) 


1 


A-17 


